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SYNOPSIS 


Among the numerical methods used for transonic flow 
computation, the Integral Equation Method (IEM) is one of the 
earliest. The non-linear integral equation for symmetric flow 
past a thin two-dimensional airfoil governed by the Transonic 
Small Perturbation Equation (TSPE) was first derived using 
Green’s theorem in 1950. The IEM has since then been used for 
computing approximate numerical solutions to many transonic 
flow problems. 

In this thesis entitled ' Numerical Solutions to 
Transonic Flows Past Airfoils and Alender Bodies by the Integral 
Equation, Method* , we are concerned with the application of the 
IEM for obtaining acceptable super-critical transonic f low 
solutions with shocks. The first problem considered is that of 
axisymmetric transonic flow past slender bodies of revolution 
with wind tunnel wall interference effects, satisfying the TSPE 
and linearized tangency boundary condition on the body surface 
and appropriate boundary conditions on the tunnel wall. Solu- 
tions were computed for general parabolic bodies of revolution 
by solving the non-linear integral equation numerically. 

The success, achieved in this attempt led to the idea 
that the IEM can be applied easily to obtain two-dimensional 
full-potential solutions for super-critical flows past arbitrary- 
airfoils satisfying the exact boundary conditions. The present 
approach to this problem combines the idea of internal singu- 
larities for incompressible flows with the notion of a field 
source distribution representing the effects of compressibility 



xvii 


to formulate a system of two non-linear integral equations for 
the velocity components * Introducing artificial viscosity and 
upwind differencing, a simple direct iteration scheme is employed 
to obtain, super-critical solutions with shocks for flows past 
arbitrary airfoils. 

Chapter I of the thesis gives a review of the various 
numerical methods used at present for transonic flow calculations 
employing the potential flow equations. The Integral Equation 
Method, in particular, is reviewed in detail and the chapter 
concludes with a discussion of the motivation and scope of the 
thesis. 

Chapter II presents the formulation of the problem for 
axisymmetric, transonic, wall-interference flow past slender 
bodies employing the small perturbation equation. The appli- 
cation of Green* s theorem to the flow domain leads to the 
derivation of an integral equation for the axial perturbation 
velocity component. The integral equation is solved numerically 
by appropriate iterative procedures. The chapter concludes with 
a discussion of the computed solutions for flows past general 
parabolic bodies of revolution in free air as well as with 
porous will tunnel wall interference. 

Chapter III presents the formulation of the problem 
for plane two-dimensional transonic flow past arbitrary airfoils, 
using the full-potential equation and exact tangency boundary 
condition. Two simultaneous non-linear integral equations are 
formulated for the two velocity components by considering the 
velocity induced at a field point by a distribution of singu- 
larities on the mean line of the airfoil together with a field 
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source distribution representing the compressibility effects. 

The integral equations are discretized on a curvilinear grid and 
solved numerically by a direct iteration scheme. Computed solu- 
tions for several test airfoils are presented in the concluding 
section of this chapter. 

Chapter IV lists the conclusions to be drawn from the 


present work. 



CHAPTER I 


INTRODUCTION 


1.1. GENERAL 

The high lift to drag ratio of wings for free -stream 
Mach numbers just below the drag -rise Mach number [1] has long 
been recognized by aerodynamicists to make transonic flight 
design conditions an attractive proposition. However/ it is 
precisely in this regime too that the flow has been recognized 
to be highly unstable. Early wind tunnel tests had already 
revealed that as the free-stream Mach number approached unity, 
the flow tended to become unstable leading to violent fluctu- 
ations in the surface press lire distribution and shock loca- 
tions [ 2 ]. This posed problems for the stability of the air- 
craft crossing the transonic regime in flight. Compression 
shocks, which are typical features of transonic flow, may be 
strong enough to cause separation of the boundary layer leading 
to loss of lift and stalling of the aircraft. The high 
sensitivity of transonic flews to small changes in the flow 
parameters makes transonic wing design difficult. Detailed 
and accurate knowledge of the flew field throughout the trans- 
onic speed range is essential to specify reliable design and 
actual flight conditions. This can be acquired by either 
wind tunnel testing of models or computational schemes based 
on mathematical modelling of the real f low. In practice a 
judicious combination of both appears to be necessary. 
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Before the advent of the modern computers , aerody- 
naroicists relied mainly on wind tunnel data for studying the 
flow fields around wing geometries for high subsonic and 
transonic speed ranges. Wind tunnel testing has become very 
expensive now. It was also recognized that wind tunnel simu- 
lation of unbounded flow over the test body was seriously 
impaired by wall interference effects at high subsonic free- 
stream Mach numbers. Thus for high subsonic and transonic 
flows, the development of analytical or numerical methods of 
solution was particularly important. Analytical results for 
transonic flows were obtained by using approximate methods to 
solve the Transonic Small Perturbation Equation (TSPE) for 
simple cases. However it is the development of faster and 
larger computers in the sixties and seventies, which made 
feasible for the first time, the computation of entire flow 
fields by numerical methods. The cost of computation too 
decreased remarkably, in inverse proportion to increasing 
computer capabilities [3], This gave the necessary impetus 
to computational transonic aerodynamics in preference to costly 
wind tunnel testing. 

The computational methods for transonic flows may be 
classified according to the choice of the governing equations 
and the method of obtaining solutions. The governing equations 
generally used are (i) the Navier-Stokes equations, (ii) the 
Euler equations and (iii) the full -potential equation or the 
TSP equation. The methods of solution may also be classified 
under each of these headings. The three approaches differ 
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basically in the idealizing assumptions made about the real 
flow which then lead to the appropriate model equations. 

In in viscid transonic aerodynamics/ the chief 
feature of the flow is the appearance of weak shocks across 
which entropy rise is negligible. Crocco-Vas zony i equation 
then states that such flews are also irrotational so that a 
velocity potential exists. Hence the equations describing 
inviscid compressible flow may be reduced to basically two 
equations/ that of mass continuity and energy conservation. 
Although such a simplification is of enormous computational 
value# the results obtained by this approach/ in the absence 
of appropriate entropy information, lead to uncertainties 
about the location and strength of calculated shocks. Calcu- 
lations by Tai [4], using the Euler equations have shown that 
some of these uncertainties may be explained on the basis of 
Euler solutions assuming appropriate entropy increments across 
the shocks. In this thesis however, we shall confine ourselves 
to potential flow solutions obtained by numerical methods-, 
keeping in mind their relevance to the Integral Equation 
Method (IEM) . 

A brief review of some of the methods that are widely 
used for computing transonic flows, is presented in this 
chapter. One of these, the integral equation method is 
employed in this thesis for computing axi symmet ric and plane 
two-dimensional flows. A review of the IEM and a discussion 
of the scope of the thesis are presented later in the 
chapter. 
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1.2. FINITE-DIFFERENCE METHODS ( FDM) 

Most of the research efforts in transonic flow compu- 
tation during the last decade has gone in to the development 
of finite-difference calculation procedures for two-dimensional 
and three-dimensional flow* past wings/ wing-body combinations 
etc. The pioneering work of Murman and Cole [5] in 1971/ 
which virtually initiated the era of modern computational 
transonics , belongs to this category. They used a type- 
dependent , non -conservative differencing scheme to solve 
successfully the non-linear, type-mixed differential equation 
that describes small perturbation transonic flow past thin, 
two-dimensional, plane airfoils. Point and line relaxation 
techniques were used to obtain iterative solutions to the 
system of non-linear algebraic equations resulting from the 
differencing scheme. Shocks, if any, emerged naturally in the 
course of the iterative procedure, as regions of continuous 
flaw with steep velocity gradients. Since then, relaxation 
schemes have been applied with great success to a wide range 
of transonic flow problems. 

The mixed differencing scheme of Murman and Cole (1971) 
suffered from the drawback, that it led to the use of an elli- 
ptic, centred difference operator at the first subsonic point 
downstream of the shock. This resulted in the shock jumps 
and shock locations being incorrectly calculated, the calcula- 
ted strength being smaller than that predicted by the conser- 
vation laws. This deficiency was rectified by Murman in a 
later paper [6] in which the shock point operator was separately 
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introduced as the sum. of the elliptic and hyperbolic operators 
and this was shown to satisfy the conservation laws. It also 
led to the correct jumps in the actual calculations . 

The general method for constructing consistent diff- 
erence schemes for hyperbolic systems of conservation laws 
had been given earlier by Lax and Wendrof f [ 7 ] . They had shown 
that the correct jumps would be calculated only if the conser- 
vation form of the equations was retained in the differencing 
scheme. The mixed differencing scheme of Merman and Cole was 
later put on a firmer footing by Jameson [8] and Albone [ 9] , 
who introduced the rotated difference scheme whiah makes the 
differencing scheme virtually Independent of the co-ordinate 
system by correctly including the domain of dependence of the 
differential equation in that of the difference equation (the 
Courant -Friedrich -Lewy condition). Thus with the use of 
conservative differencing together with type-dependent opera- 
tors and rotated differencing at supersonic points , the problem 
of computation of inviscid transonic flows may be said to have 
been solved completely in principle. Later efforts in compu- 
tational transonics have gone mostly into the development of 
algorithms that are faster and more stable and also in deve- 
loping satisfactory solution procedures for complicated three- 
dimensional geometries. 

The "speed of computation" problem has been approached 
in several ways. In the initial stages of the development of 

relaxation schemes* acceleration techniques akin to the 

2 

Aitken 6 -process, developed by Hafez and Cheng [10] were 
successful in achieving, at least in some cases, an order of 



6 


magnitude reduction in the number of iterations. The success 
of these techniques was however limited to a restricted class 
of flows. Other approaches incorporating hybrid iterative 
schemes have been tried by Jameson [ll] , giving encouraging 
results. In this method a fast Poisson solver is used after a 
specified number of relaxation sweeps to achieve faster conver- 
gence. The Poisson solver is unstable in the supersonic region 
and the relaxation sweeps are necessary to prevent instabil- 
ities from amplifying . The use of this technique is however 
limited by the strict conditions imposed on the grid spacing 
and the treatment of the boundary conditions upon which the 
efficiency of a Poisson solver critically depends. An alter- 
native approach to the speed problem has employed the 
Approximate Factorization (AF) algorithms. The idea underlying 
this approach is that, if the linear operator used to evaluate 
the correction vector in an iterative scheme closely approxi- 
mates the discretized partial differential operator vised to 
evaluate the residual, the iterative scheme will converge 
rapidly. The AF algorithms seek to achieve this by expressing 
the linear qperator as a product of factors each of which is 
easily invertible and which is a better approximation to the 
discretized partial differential operator than the unfactored 
one. A series of papers published in recent years indicate 
that this approach is superior to the line or point relaxation 
schemes [12-16]. 

A third approach engaging the attention of researchers 
has been the so called Multi Level Adaptive Technique (MLAT) 
[17-19]. In the MLAT, the discretization and solution processes 
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are Intermixed . A sequence of uniform grids or levels with 
geometrically decreasing mesh sizes, participates in the 
process. The co-operative solution process on these grids 
involves relaxation sweeps over each of these grids, coarse 
grid to fine grid interpolation of corrections and fine grid 
to coarse grid transfer of residuals. Brandt and South [20] 
have applied the MLAT to the solution of TSPE and observed a 
high rate of convergence. Jameson [21 ] and McCarthy and 
Reynher [22] have recently used the method for transonic 
potential flew calculations , the former for flow past airfoils 
and wings and the latter for flow in axisymmetric inlets. 
Preliminary results indicate a substantial increase in the 
efficiency of the calculations. Some of the earlier diffi- 
culties experienced in the treatment of lifting flows and in 
calculations on non-uniform and curvilinear meshes have been 
tackled successfully by combining the multi-grid method with 
a generalized alternating direction method [21] . 

Falling within the category of finite-difference 
methods but differing from it by virtue of employing direct 
solvers, the semi -direct method developed by Martin and Lomax 
[23-24], appears to be a promising alternative for tackling 
the speed problem. In this approach the governing differential 
equation is formulated as a Poisson equation with the non-linear 
terms acting as forcing functions to which Murman * s mixed 
difference operators are- applied. The left side of the equation 
retains its elliptic character which, in the case of super- 
critical flews, is balanced by appropriate terms added to the 
right so that at convergence, the elliptic character inducing 
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terms cancel out on either side. Jameson [25], while giving 
an excellent review of the finite-difference methods for tran- 
sonic flows, has pointed out that such an elliptic operator on 
the left is unstable for purely supersonic flows unless desym- 
metrized by an additional term. The addition of this desymm- 
etrization term however destroys the regular structure required 
by most direct solvers. This method is of particular interest 
to us because the I EM casts the governing equation in the 
Poisson form and transforms it into an integral equation. Thus 
while sharing with the semi -direct method the disadvantages 
arising out of the use of an elliptic operator, it differs 
from it in the implementation of the body surface and far-field 
boundary conditions which are now expressed as integrals over 
appropriate bounding surfaces. This implies that the boundary 
conditions in the IEM are satisfied relatively more simply and 
accurately . 

1.3. FINITE-ELEMENT METHODS (FEM) 

Subsonic flow around airfoils and wings, formulated as 
elliptic boundary value problems, have been calculated using 
finite-elements by Argyris and Mareczek [26] and Periaux [ 27 ]. 
The presence of embedded supersonic zones and shock waves 
complicates the application of finite-element methods to the 
computation of transonic flews. 

The basic idea common to most finite-element methods, 
is that the problem of finding a solution to a boundary value 
problem in partial differential equations is reformulated as 
the problem of minimisation of a certain functional which may 
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or may not be subjected to constraints. In their review of 
the finite-element methods in transonic flews , Hafez et.al. 

[28] have classified these as hyperbolic , mixed and elliptic 
methods. 

The hyperbolic method is applicable to hyperbolic 
partial differential equation systems. An approach using the 
Galerkin weighted residual technique was applied to the time- 
dependent transonic small perturbation equation by Wellford and 
Hafez [29]. Phares and Kneile [30] have applied a similar 
method to the equations of momentum conservation. Some preli- 
minary results have been obtained using these methods but the 
convergence rate is rather disappointing compared to, for 
example, the finite -difference methods. 

The mixed type methods are a finite -element analogue 
of backward differences, artificial time. and artificial visco- 
sity in the finite -difference methods. This method has been 
used by Hafez et.al. [3l] who have calculated flows over a 
parabolic airfoil using rectangular elements and line relaxa- 
tion technique. Chan and Brashears [32] have adopted a least 
squares approach to seek solutions to the TSPE by minimising a 
certain functional. This method fails to converge unless up- 
stream effects are cancelled at supersonic nodes through a 
special assembly procedure. 

Extensions of these methods to the full -potential equa- 
tion are not trivial because here the flew direction is comple- 
tely unknown. Difficulties therefore arise in the application 
of finite-element methods to super-critical flows involving 
change of strategy across the sonic line and shock wave. 
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Among the elliptic methods/ Wellford and Hafez [33] 
adopt a variational approach and solve a minimisation problem 
for full-potential flow by standard finite -element techniques. 
They used a dual iteration method in which a direct Poisson 
solver is used for the potential f unction and a gradient method 
for the velocity function. An artificial viscosity was added 
to make the scheme converge- Results were obtained for a 
parabolic arc airfoil at super-critical Mach numbers. 

Glcwinski and Pironneau [34] have adopted an optimal 
control approach and solve a minimisation problem subject to 
constraints specified by the TSP or FP equation. They expli- 
citly imposed the condition that the normal velocity down- 
stream of the shock be less than that upstream and obtained 
results for symmetric flow past a NACA 0012 airfoil. 

Eberle [35] has applied the finite-element method to 
the problem of airfoil optimisation in transonic flow. The 
addition of artificial viscosity enables the procedure to 
be converted into an analysis method. Solutions have been 
obtained by Eberle for flows past two-dimensional airfoils. 

More recently Habashi and Hafez [36] have applied the 
finite element method to solve the two-dimensional/ full- 
potential equation. Using Galerkin's weighted residual tech- 
nique/ and using artificial compressibility / they have obtained 
symmetric solutions for flows past NACA 0012 airfoil at 
subsonic and sxper sonic free -stream Mach numbers. 

The Finite-Element method has also been used to obtain 
transonic flow solutions to the full-potential equation by 
reformulating the latter as an optimal control problem. Using 
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this approach, Bristeau [37] has obtained super-critical 
solutions for symmetric and lifting flows. 

Although the finite-element method has not yet become 
competitive with finite-difference methods for transonic flow 
calculations, for problems involving geometrical complexity 
there is reason to believe that it will become so. For the 
TSP equation, using linear shape functions on triangle elements 
and line relaxation methods, Hafez et.al, [30] shewed that 
finite-elements are comparable to Wurman 1 s finite -difference 
calculations. The validity of this conclusion in the case of 
the full-potential equation however remains to be seen. A more 
promising approach appears to be the finite -volume approach 
which is in fact a combination of the finite -element and finite- 
difference approaches and to which we will now turn our atten- 
tion. 

1.4. FINITE-VOLUME METHODS ( FVM) 

A major difficulty experienced in using finite- 
difference methods is that associated with the implementation 
of boundary conditions for different body geometries. For a 
Cartesian mesh, special finite-difference formulas are needed 
near the boundary leading to loss of accuracy. An alternative 
approach is to employ mappings to obtain geometrically simple 
boundaries on the computational domain. This may sometimes 
be a formidable task, especially for realistic three-dimensional 
geometries. A third approach, circumventing these difficulties, 
is provided by the finite-volume technique. 
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Finite- volume methods use general, non- orthogonal 
co-ordinates and consider the governing differential equations 
as balances of mass, momentum and energy fluxes for each 
finite- volume defined by the intersection of co-ordinate 
surfaces. This method was first introduced by MacCormack 
and paullay [ 38] who studied the discontinuous or weak solu- 
tions of the time-dependent two-dimensional Euler equations 
for supersonic flow past a wedge, a double symmetric wedge 
and a sphere. Rizzi [39], Caughey and Jameson [40-42] have 
all applied the finite- volume method to representative tran- 
sonic flow problems. Caughey and Jameson use local bilinear 
mappings to obtain curvilinear cells on the physical plane 
from regular cubes on the computational plane. A similar 
bilinear variation, is assumed also for the velocity poten- 
tial. They then use the usual finite-difference techniques 
to cast the quasi-linear equations in difference form, while 
ensuring that the flux balances across element interfaces 
satisfy the mass conservation, law. 

The finite-volume method combines the advantages 
of finite-elements for handling complicated geometries 
with the simplicity of finite-differencing. Jameson and 
Caughey [42] have reviewed the progress in the applica- 
tion of finite-volume methods to calculation of tran- 
sonic full-potential flows past general three-dimensional 
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wing-body combinations. They observe that though the solid 
surface boundary conditions are easily satisfied, satisfaction 
of far -field boundary condition suffers from the need to 
truncate the boundary at a finite distance from the body. 

This is a disadvantage that finite-volume methods still share 
with finite-difference techniques. Nevertheless, for many 
three-dimensional configurations of aerodynamic interest, grid 
construction is considerably simplified by the use of the 
finite-volume technique which requires co-ordinate transfor- 
mations to be defined only locally. The accuracy in imple- 
menting the surface boundary condition is thereby improved 
without making the algorithm, unreasonably complex. 

1.5. INTEGRAL EQUATION METHODS (IEM) 

In this thesis we are concerned with the use of IEM 
for the computation of transonic flows. A comparatively 
detailed review of this method is therefore given in this 
section . 

The integral equation method is one of the earliest 
methods used to obtain solutions for transonic flews past air- 
foils, wings and bodies of revolution. In the IEM we deal 
with the flow directly on the physical plane. The potential 
flow equations are transformed into integral equations by the 
use of Green's theorem. This gives the perturbation potential 
at a general point in the flow field in terms of line or 
surface integrals over the flow boundaries (solid boundaries 
and flow discontinuities treated as internal boundaries) and a 
volume integral of all the non-linear terms limped together and 
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acting as a field source distribution. The integral equation 
is then discretized by dividing the flow field into panels or 
cells and the resulting system of non-linear algebraic equations 
is solved iteratively. 

Most of the approaches to the transonic flow problem 
by the IEM to-date have used the transonic small perturbation 
equation , together with the linearized, mean-line/mean-surface 
approximation for satisfying the tangency boundary condition. 
Among the earliest such attempts, we must mention three related 
studies by Oswatitsch, Gullstrand and, Spreiter and Alksne 
(hereafter referred to as OGSA approach). Oswatitsch [43] 
deduced a non-linear integral equation for a profile at zero 
lift, resulting in a line integral around the airfoil and a 
two-dimensional space integral. He treated shock waves as 
internal boundaries, approximately normal to the airfoil surface 
and across which the Rankine-Hugoniot jump conditions are satis- 
fied. Introducing an assumed variation of the axial velocity 
component in the transverse direction, he reduced the two- 
dimensional problem to one with a single independent variable. 
This problem was then solved by quadrature. This work was 
later extended by Gullstrand [44,- 48] who also derived the 
integral equation for symmetric flow past three-dimensional 
wings and examined planar two-dimensional lifting flows. 

Spreiter and Alksne [49] solved the integral equation for 
sub-critical and continuous super-critical flews by iteration. 
For flows with shocks, the shock location on the airfoil was 
prescribed and the free -stream Mach number determined. In all 
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these cases, the method was restricted to two-dimensional 
symmetric flows past sharp nosed airfoils. 

The method of Oswatitsch was later extended by Keune 
and Oswatitsch [50 ] to include the ax;i sy mme t r i c case, Heaslet 
and Spreiter [51] applied the I EM to transonic flow around 
slender wings and bodies of revolution. Special attention was 
given to the conditions resulting from the presence of shock 
waves and to the reduction of relations to special forms 
necessary for the discussion of flows at Mach one. Results 
were obtained for cone cylinders, wings and wing -body combina- 
tions at unit free stream Mach number and compared with experi- 
mental results. Klunker [52] deduced the integral equation for 
the velocity potential for transonic flew past thin lifting 
wings in subsonic f ree-streams, treating shock discontinuities 
appropriately. He also obtained the analytical expression for 
the far-field. 

Np'rstrud [53-56] extended the OGSA approach to lifting 
and non-lifting flows with shock discontinuities. The diff- 
erent f uncticnal approximations that can be made regarding the 
transverse variation of the asq.al velocity component were 
examined [53], Applying the method of parametric differentia- 
tion, a Cauchy initial value problem was setup and solved by 
the numerical techniques of ordinary differential equations. 
This method however works only for shock free flows. For flews 
with shock discontinuities, the method of steepest descent was 
used. The initial velocity distribution in this case was 
assumed to have a discontinuity at an appropriate station. The 
solutions obtained were found to compare favourably with 
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finite “diffe r enc e and experimental results. These techniques 
were extended to cover three-dimensional finite lifting wings 
[55-56], Nixon and Hancock [57] derived two simultaneous 
integral equations for two-dimensional, planar lifting flows. 
Computed results showed close agreement with exact results [ 112] 
obtained by Sells. The method was applied by Nixon [58] to an 
airfoil oscillating at low frequency in a high Mach number 
subsonic flow. 

Niyogi and his associates used the basic method of 
Oswatitsch to obtain approximate analytic solutions which are 
then used as input in an iterative solution scheme. These 
studies have been reported by Niyogi [59] in a review article. 
Niyogi [60] also sought to apply these techniques to three- 
dimensional finite symmetric wings at zero incidence and sub- 
critical speeds. 

All the studies mentioned so far have dealt with the 
problem for free air flow only, Kraft and Lo [61] used the 
IEM for determining the effect of a perforated wind tunnel wall 
on the flow past a two-dimensional, non-lifting, plane airfoil 
at transonic speeds. They demonstrated that the correct shock 
location of free air flow can be obtained by the proper selec- 
tion of porosity. Earlier Mokry [62] applied the IEM based on 
linear theory to compute subsonic flows past lifting two- 
dimensional airfoils in the presence of ventilated walls, 
assuming a general form of the classical wall boundary condition 
to simulate open jets, porous walls and slotted walls. 

These early investigations using the IEM were severly 
restricted by computer capabilities available during the 
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sixties end early seventies. This limitation forced the 
investigators not to deal directly with the two-dimensional 
integral equation but to reduce it to one dimension/ lasing 
arbitrary functional approximations for the transverse variation 
of the axial velocity component. This assumption is a feature 
common to all the methods discussed above. It is true that such 
a functional approximation/ while drastically reducing the 
computational effort/ does not seriously affect the accuracy of 
the final results at least for shock-free flows. For super- 
critical flows with shock however/ the monotonic behaviour of 
the velocity inplied by such functional approximations/ breaks 
down in the vicinity of the shock, thus invalidating the whole 
procedure. Modifications to rectify these errors were suggested 
by Nixon [63] who introduced the ‘extended integral equation 
method' wherein the flow field is divided into a finite number 
of streamwise strips and the transverse variation of the pertur- 
bation velocity across each of these strips is approximated by 
interpolation functions in terms of the values at the strip 
edges. Nixon [64-67] showed that this modification gave consi- 
derably improved results. 

The first attempt to find exact sub-critical numerical 
solutions to the two-dimensional integral equation without any 
approximating assumptions, was made by Ogana and Spreiter 
[68-69], In this approach, the flow field was divided into 
rectangular panels and the discretized system of non-linear 
algebraic equations was solved by Jacobi iteration [69] , Results 
were obtained for sub-critical flows past symmetric airfoils. 
Arora and Agarwal [70] used the field panel technique for the 
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problem of sub-critical axisymmetric flow past slender bodies 
of revolution while Arora [7l] treated transonic flow past 
lifting airfoils including wind tunnel wall interference effects. 

Lobo and Sachdev [72] used higher order functional 
approximations for the accurate numerical evaluation of the 
singular field integral. They solved the problem of disconti- 
nuous super-critical flow past symmetric profiles at zero angle 
of attack, by the quadratic iteration scheme. 

Pier® and Slooff [73] obtained solutions for super- 
critical transonic flow past a symmetric parabolic airfoil by 
the explicit introduction of an artificial viscosity term into 
the TSPE. The field source term was subjected to conservative 
differencing as in Jameson's scheme. A field panel approach 
with a point source approximation for the square field panels, 
was used to discretize the integrals and the resulting system 
of equations was solved iteratively by quasi-Newton scheme. 
Calculated results for a parabolic airfoil show good agreement 
with the results of finite-difference calculations. Similar 
ideas have been pursued by Lobo and Sachdev [74] who used 
artificial viscosity to obtain smeared shock solutions which 
were then used as the initial input to a quadratic shock-fitting 
scheme, yielding sharp shocks satisfying the Rankine-Hugoniot 
jump relations. An approach similar to that of Piere and 
Slooff was used by Ravichandran , Arora and Singh [75] for axi- 
symmetric flow past slender bodies in the presence of perforated 
wind tunnel walls. Artificial time-dependent terms modeled on 
Jameson's scheme [8] were added to impart stability to the 
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iterative scheme. This effort will be described in detail in 
the following chapter. 

All the methods described so far employed the small 
perturbation equation together with the linearized/ mean-line 
approximate boundary condition as the governing relations for 
obtaining the integral equations. An important advantage of 
the small perturbation theory is that it permits the application 
of the tangency condition on the body axis for two-dimensional 
and axisymmetric flows and on the plane of the wing for three- 
dimensional flows. In a finite-difference calculation proce- 
dure, such an approximation facilitates the use of a uniform 
Cartesian grid so that the need for special formulas or inter- 
polations to satisfy the boundary condition on the body surface 
when the latter does not lie on a co-ordinate line or surface, 
is obviated. But the mean surface approximation is obviously 
not valid near the nose region of an airfoil. Moreover, the 
fact that most of the flow problems, even in inviscid aerodyna- 
mics, are of the singular perturbation kind, implies that the 
TSPE solutions are not uniformly valid throughout the region of 
interest. Further, since the small perturbation equations are 
dependent upon the choice of the perturbation parameter, their 
applicability is confined to a restricted class of flows and 
body geometries. Small perturbation theory is often "tuned 
up" by scaling the variables in the problems by powers of the 
free-stream Mach number to obtain more accurate approximations 
to solutions of the full-potential equation [76>]* For super- 
critical flows however, the location and strength of the shock 
may be strongly dependent upon the values of the similarity 
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exponents. These exponents are not known apriori for specific 
cases and in a sense they may even be said to depend upon the 
solutions to be obtained [77] » Hence attempts to tackle more 
general flow problems by the small perturbation approach are 
beset with many difficulties and restrictive assumptions. 

Finally , recent bench mark data obtained by Hinson and Budges 
[78] have led to the conclusion that “full-potential solutions 
are uniformly in better agreement with experiment than small 
disturbance results 11 . The development of larger and faster 
computers has made it possible to dispense with the small 
disturbance theory altogether in favour of the full-potential 
equation and exact boundary conditions. This approach has also 
made it possible to develop computer codes capable of handling 
a wide class of flows so long as they satisfy the weak shock 
assumption for i sen tropic , irrotational flows. 

Exact integral equation formulation for sub-critical 
flows using the full -potential equation has been studied by Luu 
and Goulrny [79] and by Strieker [80] . Combining the classical 
surface singularity approach with the field panel method, 
Strieker has obtained solutions for sub-critical flows past 
lifting airfoils.. Comparison with Sells* [112] results show good 
agreement. These methods however do not cover supercritical 
flows with shocks. In the present work, an exact integral 
equation formulation for transonic flow using the full -potential 
equation has been employed for the computation of super-critical 
solutions with shocks for flaws past two-dimensional airfoils. 
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1.6. MOTIVATION FOR THE PRESENT WORK AND SCOPE OF THE THESIS 

In the light of the discussion above, it appears that 
the integral equation method could be a promising alternative 
for the computation of at least some classes of transonic flow 
problems, particularly those involving complex geometries. In 
viewing the I EM as a possible alternative to the finite -difference , 
finite -volume and finite-element methods, one has two prime 
consideration in mind: (i) easy and accurate implementation of 
the boundary conditions and (ii) the possibility of obtaining 
reasonably accurate solutions on relatively coarse grids in 
fewer number of iterations. 

It has already been observed in the discussion of the 
finite-difference and finite -volume methods, that the main problems 
encountered therein concern the implementation of the boundary 
conditions. In flow problems the boundary conditions are usually 
of three types. One is the body surface boundary condition.. 

This is usually tackled by globally mapping the physical domain 
on to a rectangular domain as in finite -difference methods or 
by using locally defined mappings as in finite-volume methods. 
Although this results in boundaries conforming to co-ordinate 
lines or surfaces in the computational domain, making it possible 
to use simple difference formulas for the derivatives, the 
necessity of using one sided (backward or forward) differences 
may lead to loss of accuracy in the calculation of normal deri- 
vatives on the boundary points. To overcome this difficulty, 
extremely close spacing of grid lines in the vicinity of body 
boundaries becomes essential. The second problem is related to 
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the far -field boundary condition. This is usually satisfied by 
prescribing an analytic far -field approximation at the outer 
boundary of the grid on the physical plane, or the entire 
physical domain is mapped on to a finite computational domain in 
which case the vanishing far-field conditions are directly 
satisfied. In the finite-volume approach where only local 
mappings are employed, the far-field conditions are assumed to 
prevail at a finite distance away from the body, to which the 
computational domain is then arbitrarily restricted. A third 
type of boundary condition sometimes arises across internal fluid 
boundaries such as wakes or trailing vortex sheets but we will 
not be concerned with them here. 

The I EM, in this respect, offers an attractive alter- 
native. Here the boundary conditions are implicitly satisfied 
by the integral equations in the form of integrals over bounding 
surfaces. Since integration is an error smoothening process 
and therefore inherently more accurate than that of finding 
derivatives by talcing differences, the body surface boundary 
condition can be accurately satisfied more or less independently 
of the grid used for the overall computation. Also since the 
far-field boundary condition is implicitly contained in the 
integral equations, the computational domain can be justifiably 
limited to a finite distance around the body. Though co-ordinate 
transformations may be employed, as done in the present work, 
they merely serve to set-up a convenient computational grid and 
are used only for the purpose of evaluating the derivatives of 
flow quantities. Thus in terms of accuracy, flexibility and 



ease of implementation the method appears to be of promise 
atleast for two-dimensional flows. 

Some of the problems associated with the use of the 
IEM may also be recognized at the outset. One of these is the 
need to compute and store a large number of influence co- 
efficients. The consequent severe restriction on the number of 
grid points , though not a problem for sub-critical calculations, 
could prove unsatisf actory for the computation of super -critical 
flows where sharp shocks may have to be captured. Secondly, 
instability in the numerical scheme for strongly super-critical 
flews may be expected by the use of a Poisson type of iteration 
scheme which suggests itself naturally for the present formu- 
lation. It should be noted that such an iteration scheme has 
failed to produce convergent results even for slightly super- 
critical flows in finite-dif f erence calculations [23] . 

In the present work the IEM has been used to find TSP 
and FP solutions to axisymmetric flow past slender bodies of 
revolution and plane flow past two-dimensional airfoils respec- 
tively. In this approach, an artificial viscosity term which 
enables shock capture is added explicitly to the governing 
differential equation. Non-linear integral equations for the 
perturbation velocity components are then formulated. The 
integral equations are discretized on a rectangular or a curvi- 
linear grid. The discretized system of equations is solved 
numerically \osing an appropriate iteration scheme. 

In dealing with the IEM, the present thesis is concerned 
with the application of the method for the computation of prac- 
tical flows. Theoretical questions regarding the convergence 
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and stability of the numerical schemes and the satisfaction of 
conservation laws across shocks will not be dealt with. Further/ 
though the thesis gives a good idea of the computational effort 
involved/ a comparative evaluation of the IEM with other methods, 
particularly finite -difference methods, has not been provided 
for obvious reasons. Such an evaluation is of course necessary 
in the long run but of questionable value at present because, 
to quote Golberg [81 ] , "where comparisons are made ’first-cut ' 
programs for integral equations are put up against 'state of 
the art* codes for partial differential equations. Even so the 
integral equation may win out " . In this sense , the present work 
is to be regarded still as an exploratory first-cut attempt at 
using the IEM for transonic flow calculations. 

In what follows, in Chapter IT we present the formulation 
as well as the results of the problem for axisymmetric flow past 
slender bodies of revolution, using the transonic small pertur- 
bation equation and the linearized body boundary condition. 
Chapter III presents the full -potential formulation and the 
solutions obtained for plane flow past two-dimensional airfoils. 
Chapter IV, the last, lists the conclusions to be drawn from the 


present work 



chapter II 


AXISYMMETRIC TRANSONIC FLOW PAST SLENDER 
BODIES -SMALL PERTURBATION SOLUTIONS 

2.1. INTRODUCTION 

In. -this chapter, the axisymmetric transonic small 
perturbation equation, together with the linearized boundary 
condition on the body surface and appropriate outer boundary 
conditions for wind tunnel or free air flow, is converted into 
an integral equation for the axial perturbation velocity conpo- 
nent by the application of Green's theorem. The integral equa- 
tion is then discretized into a system of non-linear algebraic 
equations by dividing the flow field in the meridian plane into 
rectangular panels wherein the flew variables may be assumed 
constant. Using a quasi-Newton or a direct iteration scheme, 
the system of non-linear algebraic equations is solved to obtain 
surface pressure distributions for slender, parabolic bodies of 
revolution for free air and perforated wall cases. Section 2.5 
of this chapter discusses the results of the computations and 
the performance of the iteration schemes for the axisymmetric 
case. The chapter concludes with a brief mention of some 
TSP results obtained for plane symmetric flaw past parabolic 
arc airfoils. 

Before presenting the formulation of the problem, a 
few remarks of an introductory nature will be made on the 
question of wall interference in transonic flow testing. 
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2.1*1. Wall Interference Effects 

The object of wind tunnel testing of aerodynamic 
configurations is to acquire data that would predict the perfor- 
mance of test bodies (wings , fuselages# wing -body combinations 
etc.) under actual flight conditions in an infinite# interfer- 
ence free domain. It has however been recognized since long 
that# even for low speed conditions# the finite test section 
size of wind tunnels produces effects on the flow around the 
test object which result in the measured pressure distribution 
being different from, what it would be in free flight. Along 
with these effects, another typical feature of transonic flow 
also comes into play — that the disturbances in the transverse 
direction decay much more slowly and hence extend much farther 
than for low subsonic flows# thereby demanding a larger tunnel 
section. Thus for super-critical flows with embedded super- 
sonic pockets# the effect of tunnel walls becomes more prono- 
unced and must be accounted for. Hence arises the need for 
precise quantitative and qualitative evaluation of wall inter- 
ference in wind tunnel testing in general and transonic wind 
tunnel testing in particular. 

There are basically two approaches to the problem of 
wall interference in transonic wind tunnel testing. One is the 
adaptive wall concept [82]# the basic idea of which is as 
follows: In any wind tunnel test of a model configuration # it 

is possible to measure# jit a convenient control surface away 

' .S' . . 

from the model and say near the walls# two perturbation flaw 
quantities# such as flew deflection and static pressure or 
axial and transverse velocity components . Using one of the 
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measured quantities as the boundary condition on the chosen 
control surface, the exterior flow is solved for by satisfying 
the governing equations. This yields a prediction of the other 
flow property. The predicted value is compared with the 
measured value to assess their compatibility. The wall condi- 
tions are then changed iteratively until the computed and 
measured values agree, implying interference free flow in the 
tunnel. The second approach is the correctable interference 
concept [83 ] in which the distribution of a single experimentally 
measured flow quantity along a conveniently chosen control 
surface far from the model is used as the outer boundary condi- 
tion in a computational procedure to obtain the flew field around 
a model. The comparison of these results with the measured 
values given an assessment of non-potential effects, while 
comparison with free air computations indicates the presence of 
wall interference. The correctable interference concept hence 
obviates the need to model the tunnel wall boundary condition. 
However there is no way of applying precise quantitative correc- 
tions to tunnel data for obtaining interference free results. 

Calculations to assess wind tunnel wall interference 
have been carried out by Bailey [84] and Sedin [85] for axi- 
symmetric flows and Murman [86] and Kacprzynski [87] for two- 
dimensional flews. Bailey, Sedin and Murman have used the TSPE 
formulation together with the classical homogeneous wall boun- 
dary condition. Kaeprzynski, in addition to the TSPE formula- 
tion, has also used the full-potential equation. More recently 
Karlsson and Sedin [88] have presented a method for calculating 
transonic wind tunnel wall interference in an axisymmetric 
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slotted test section. They discuss the problem of designing 
optimum. slot shapes for obtaining interference free conditions . 
Kraft and Lo [6l] have used the IEH for studying blockage 
effects in a perforated wall transonic wind tunnel. In their 
calculation however/ an approximating assumption regarding the 
transverse decay behaviour of the axial velocity component is 
made/ which though perhaps reliable for subsonic shock-free 
flews, becomes totally unreliable for super-critical flows with 
shocks* This points to the need for considering the entire 
flew field and the full form of the governing integral equation 
in order to obtain reliable results. 

In modelling the presence of porous wind tunnel walls , 
the wall boundary condition assumes crucial inportance. There 
has been some disagreement among researchers regarding the 
ability of the classical homogeneous wall boundary condition 
to simulate the actual p rope; -ties of perforated and slotted 
walls. The advantage of the correctable interference concept 
is / as already mentioned, that it does away with the need for 
such models . Stahara and Spreiter [89] , based on the correct- 
able interference concept, have carried out extensive calcula- 
tion for axisymmetric flow past slender bodies to conclude that 
"such a procedure can be useful as a practical wall correction 
method for moderately variable or fixed geometry transonic wind 
tunnels where alteration of the tunnel wall is limited or 
impassible" . In the absence of wall data however, a wall boun- 
dary condition such as the one mentioned above has to be used. 
For the perforated wall atleast, Jacocks [90 ] found that the 
classical homogeneous model is a reasonable representation. The 
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difficulty here however is that the porosity factor is not only 
a function of the flow parameters and wall geometry# but also a 
function of plenum and test section pressures# involving non- 
linear effects. The assumption that the porosity factor is 
constant along the wall# which has been made in the present 
calculations# is hence non-representative but makes the compu- 
tation simpler and serves illustrative purposes only. The 
replacement of the constant by a slowly varying function of 
distance along the wall is a sittple matter but to deal with a 
porosity factor which is a function of wall pressures may involve 
slightly different methods for the accurate evaluation of the 
wall integral. 


2.2. PROBLEM FORMULATION 


2.2.1. The Full— Potential Equation 


We consider axisymmet ric transonic flow past a slender 
body of revolution in cylindrical co-ordinates (x^ # r^# 0^) with 
the origin at the body nose and the x-axis parallel to the free- 
stream. 


The partial differential equation for the velocity 
potential $ of a compressible# inviscid# non -heat conducting 
flow of a perfect gas with only weak shock discontinuities can 
be derived from the conservation laws as 



2 . .2 
L ) L x + (c 
1 11 


2§ „ § _ § _ 
X 1 


A ) § 

r l 


r l r l 




( 2 . 2 . 1 ) 
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where c is the sonic velocity defined by the energy equation 

C 2 - K- ---5-^ ts 2 + S 2 - 1) (2.2.2) 

1 1 

In Eqs. (2.2.1) and (2.2.2)/ all velocities have been 
made dimensionless by the free-stream velocity. The length 
scales are all referred to the body length. denotes the 

free-stream Mach number and y is the ratio of the specific heats. 

The boundary condition to be applied on the body surface 
is that the flow velocity vector must be aligned with the local 
tangent to the body surface. If the body surface is prescribed 
as 

r^ = x R(x^) (2.2.3) 

where x denotes the maximum thickness of the body/ the tangency 
condition can be written as 
S r 

= r on r 1 = TR(x.) (2.2.4) 

a Y 11 

X 1 X 1 

The outer boundary condition is to be prescribed at 
the far-field for free air flow and, on the tunnel wall and far 
rep stream for well interference flow. 

The free air far-field boundary condition can be written 
as 

2 2 

"I, •* 0 as (xf + r. ) •* a> (2.2.5a) 

X 1 r l ii 

For wall interference flow, the classical homogeneous 
wall boundary condition for perforated tunnel walls may be 
written as {91] 
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5 r + P(X 1 ) ( Q ”1) 
x. 1 ± -^ 2 . 


on 


= r 


lw 


(2.2.5b) 


where r^ w represents the wall radius and p(x^) is the porosity 
parameter which relates the pressure difference across the wall 
to the flow through the wall. pte^) in general must be deter- 
mined experimentally for a particular wall geometry and free- 
stream. Mach number. In addition to Eq. (2.2.5b)# it is assumed 
that the axial velocity component vanishes far upstream and 
downstream 

§ v 1 / §_ 0 as x, *♦ + cd 

•*1 1 


me local pressure 




c = 
p Y M 


Y 

* [$1 — — — ^ — — M ($ 2 + — 1 }}*^ ^ — 1 ] ( 2 . 2 . 6 ) 

M X 1 r l 

OD 


2.2.2. The Small Perturbation Equation 

For flew past slender bodies at free-stream Mach 
numbers close to unity# the transonic small perturbation equa- 
tion may be obtained from, Eq. (2.2.1) by the method of matched 
asymptotic expansions# using an expansion parameter e which is 

a measure of the magnitude of the perturbations [92] . 

o 

Assuming that (1 - tot ) ~ 0(e)# it can be shown that 

[92] e *■» 0 (t 2 ) for axisymmetric transonic flow past a body of 
thickness Denoting the perturbation potential by 0^# the 
TSPE in physical variables can be obtained as 

[1 - >4 - (Y + 0 1X]Xi +4 = P 


(2.2.7) 
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The first order slender body boundary condition consis- 
tent with the order of approximation of Eq. (2.2.7) is 


2- dR s'Cx-) 

(r Ar, > - T R It - — ar- 

r^ **• 0 1 1 


( 2 . 2 . 8 ) 


2 —2 

where s(x^) ~ ti % r (x^) is the body cross-sectional area 
and prime denotes derivative with respect to . 

The outer boundary condition Eqs. (2.2.5) now become/ 


free air flow: 0 lx , 0 lr 


0 


as 


Cxi 


+ r?) 


1/2 


oo 

(2.2.9a) 


wall interference flow: 0 lr + p( x 2^ 0 ix 


0 on r l = r lw 

(2.2.9b) 


0 i ri ' 0 lx 1 


0 as x^ ** +oo 


Equation (2,2.7)/ toy virtue of its essential non-linearity 
and mixed character, correctly describes in viscid isentropic flow 
throughout the transonic regime. Embedded supersonic pockets 
and weak shock- discontinuities across which the entropy rise is 
negligible, are typical of super-critical transonic flews. These 
are contained in the admissible class of solutions to the non- 
linear equation satisfying the following principle of weak 
solutions 


/ 

Q 


[W i (1 

1 


)0 

OD^lX, 


Y - 4 ' - 1 M 2 rt2 


= 0 


‘I r 


( 2 . 2 . 10 ) 


for arbitrary test functions W (x^ , r^ ) . ■ q denotes the flow 

region exterior to the body. Equation (2.2.10) can be used to 
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show that the following jump condition is satisfied across a 
discontinuity [93] 


[Cl - M 2 m ) < V - \ ftf + 1)M 2 ce <U 2 >] <Uj> + < Vl > 2 


= 0 


where, u 1 = 0^.. , = 0 


!yi 


( 2 . 2 . 11 ) 


and < > denotes a jump in the quantity across the discontinuity , 
Finally the first order approximation to the pressure 
co-efficient may be obtained from Eq. (2,2.6) as 

2 


- 20 - (t 

^lx 1 dx^' 


( 2 . 2 . 12 ) 


2.2.3. Reduced Variables 

For convenience we now apply the following transformations 
to the physical variables in Eq, (2.2.7) 


x = x. 


r = pr.. 


* - hh 


(2.2.13) 


Here p 2 = 1 - M 2 ^ and K = (y + DM^, 
Eqs . (2.2.7) and (2.2.9) become 


With these transf ormatipns 


0 +0 + — 0 = 00 

^xx *rr r r x : 


xx 


lim ; (r0 ) 
r -0 r 


K_ s f (x) 

p2 2 it 


free air flows 0 ,0 

JC 


0 


, 2 ■ 2. 

as (x + r ) 


1/2 


(2.2.14) 

(2.2.15) 

<3D (2.2.16a) 


wall interference flow: 0 r + 0 X = 0 at r = pr lw = r w 

(2.2.16b) 


0 / 0 + 0 as x •* + cd 

.A. Xr 
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2.2.4. Super- critical Flow and. Artificial Viscosity 

It must be noted that Eq. (2.2.14) admits solutions with 
expansion, discontinuities which are physically inadmissible . For 
uniqueness , directional property corresponding to "entropy inequ- 
ality" must then be restored in numerical schemes.* For super- 
critical flows the directional bias can be introduced by using 
upwind differencing to represent the derivatives along the stream- 
wise direction at the supersonic points. The assumption of small 
perturbation, theory implies that the velocity vector deviates 
little from. the x~axis in the entire field. The characteristics 
in the supersonic region are then locally symmetric about the 
x-axis. Hence the use of upwind differencing in the x-direction 
ensures the correct domain of dependence at the supersonic points. 
Considering the two-dimensional TSPE in conservation form we can 
write 

f x + ay = o with f - n x - \ , a = S» y 

The error in evaluating f one mesh width Ax backward is - Axf . 

A XX 

Thus upwind differencing introduces a discretization error 
- AX -A [U - = - &x^.[(l-u)u x ] 

with u = 0 , which can be regarded as an artificial viscosity in 
the supersonic region. 

Turning this idea around, an artificial viscosity term 
of the form. 


* Inclusion, of artificial viscosity alone does not appear to 
ensure uniqueness of cortputed solutions. Multiple solutions 
have been obtained for plane flows even with the inclusion of 
artificial viscosity, see [ill] . 



35 


e [*(u)u x ] 

where e > 0 and is O(Ax)# and f (u) [S (u - 1)] > 0 and is of 
0(1) can now be explicitly added to the TSPE, representing the 
relative retardation of flow in the x -direction at the super- 
sonic points.. Equation (2.2 a 14), with the addition of artificial 
viscosity can then be written in the modified form as 

0XX ^rr + 7fr = 5<x, r) (2.2.17) 

with the non-linear source term 

g(x/ r) » uu x [f (u)u x ] (2.2.18) 

Solution of the modified Eq. (2.2.17) replace the shock 
discontinuities in the solutions to Eq. (2.2.14) by continuous 
narrow regions with steep velocity gradients. 

The choice of a suitable artificial viscosity function 
f (u) in Eq. (2.2.18)/ depends on essentially two considerations: 

(i) that the captured shock be steep given the constraints on 
mesh size and the transition through the sonic point be 

smooth and (ii) that the iteration scheme be stable and the 
convergence rate reasonable. 

For small perturbation flows past symmetric airfoils / 
several artificial viscosity functions and their asymptotic 
shock layer characteristics have been discussed by Piers and Slooff 
[ 73 ]. Among the many artificial viscosity functions tried in 
the course of the present work/ only f (u) = u - 1 was found 
suitable for axi symmetric flow computations. 
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Equations (2.2,17) and (2.2,18) together with the 
boundary conditions Eq§. (2,2.15) and (2.2.1&) completely 
specify the problem for axisymmetric transonic flow past 
slender bodies. An application of Green's theorem to the 
flow domain yields the integral equation for the perturbation 
velocity potential. 

2.3. THE INTEGRAL EQUATION 

Figure 2.1 presents a schematic representation of 
the problem specified in the previous section together with 
the domain of application of Green's theorem. 

Green's theorem states that if 0 and ¥ are two 
arbitrary functions that are differentiable twice within the 
volume Q bounded by a surf ace S then, 

/ ( 0 - 0 v 2 y )dQ » / (0 || - ® ||)ds (2.3.1) 

Q 3 3n dn 

2 

where n denotes the normal drawn into q and V denotes the 
Laplacian operator which, for the present case, in terms of 
the cylindrical co-ordinates is 

V 2 - 3 2 + a * + 1 3 + 1 $ 

For the present problem 0 is chosen as the function 
satisfying Eqs. (2.2.17), (2.2.15) and (2.2.16) while ¥ is 
assumed to be the fundamental solution of the Laplace equation 
given by 

* = l/[(x -» x' ) 2 + r* 2 + r 2 - 2rr* cos(© - ©«)] 1/2 


(2.3.2) 
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In Eqs . (2,3.2), (x, r, ©) denote the field point and (x 1 , r', 9 f ) 
denote the running variables. Note that ^ is singular when 
(x 1 , r'/ 9 1 ) coincides with (x, r, 9). 

With the above choice of 0 and W , Eq. (2.3.1) holds in 
the region Q of Figure 2.1 bounded by the surface S which 
includes the body surface S fi , an infinitesimal sphere S p surr- 
ounding the field point P(x, r, 9) and the surface S for free 
air flow or the wall surface S Tt for flow confined by tunnel 

Vv 

walls. Since shock discontinuities in the flow are assumed to 
have been replaced by continuous regions with steep velocity 
gradients, the region as defined above satisfies the conditions 
of Green's theorem. Further the flow being axisymmetric , 
without loss of generality, we may take 9 = 0 in Eq, (2.3.2). 

The various integrals occurring in Eq. (2,3.1) may be evaluated 
as follows. 

The Volume Integral: 

Since v 2 ^ 1 = 0 and ^0 = g in the region under consi- 
deration we have from Eq. (2.3.1) 

I « / (S v 2 0 - 0 v 2 ^ ) dQ 
y Q 

= /// f (x, rjr x * , r 1 , ©*) g(x", r*) r* d9 r dr 1 dx 1 

0 

(2.3.3) 

The integral in Eq. (2.3.3) being convergent, the integral 
over the spherical cavity surrounding the point P makes negli- 
gible contribution in the limit. The region q may therefore 
be taken to be the entire space. 
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The Surface Integrals s 
(±) The Integral Over Sp 

The integral over the point cavity S p may be evaluated 
as follows. The point P is surrounded by a sphere of radius £ 
centred at P, Letting e -» 0 


/ 

S P 

E •* 0 


(0 


B 1 

an 


* tS> ds 


- Lt (0(x, r) Sf |^) 4ne 2 

£ «*0 

= Lt^ [~ 0(x, r) ^ |f~] 4ite 

£ ** 0 4 7C£ 


0(x, r) 


(2.3.4) 


(ii) The Integral Over 

The integral over the body represented by the slit S g 
can be evaluated as follows 


[ S = / 

B S 


a^ 

an 


(0 377 - ® §f|) ds 


B 


When the body is sufficiently slender, the following 
approximation may be made 


lira. 


3n 


3 

3r ‘ 


Therefore 


1 2% ad 

/ / (0 ? §fr> r* d©' dx 1 


B O O 

where 1 is the location of the body sting. For bodies without 
sting 1 =» 1« In the limit r* •* 0, 0 is of the form 
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0 ® s'(x’) In r‘ + g{x‘) 


where g(x*) Is a constant of integration (see [92]), 
A1S0 from Eq. (2.3„2) 


lim 

r '•* 0 r 


r cost© - 0* ) 

l yy m i muim i m n kwiw i ii i i ay i ^ 

4tc[(x - x’) Z + r ] 3/Z 


Hence in the limit r' - 0# r" 0 ^ r , 0 and the integral 


becomes 


1 275 

« - / / 


l q - - - ¥ r* d© 1 dx* 

S B o c 9r 

r’ ■* 0 


The tangenqy boundary condition (2.2*15) may now be used to 
obtain 


'B 


1 23t 

/ / JL— gi(x') d©> dx 1 

o o 2ltp 


K 


4l5p o [(x * x*j + r* J 


s‘(x') 

.\i HOT 


dx 1 


(2.3.5) 


The integral in Eq. (2.3.5) may be recognized as giving the 
potential for linear subsonic flow past a slender body. 


(a) 


The integral Over the Outer Boundary 

For the free air case, the outer boundary condition is to 
be applied at infinity. Hence the integral over the 
surface at infinity is given by 

M 1 ! 


CD 


s 

CD 


W fs - * i£> dS 


Lit 


GO 

/ f (0 f r , - * 0 r i ) r 1 d9* dx 1 

-00 o 
r* -» oo 


(2.3.6) 
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under the assumption that 0 ^ ■i— , z > 0 and R -» od. 

R e 

(b) For the wall interference case, the integral over the 
cylindrical wall surface has to be evaluated. 


w 


J7 (0 f , - * 0 ,)r d©» dx 1 at r 1 

S x ' r w 

w 


w 


Applying the wall boundary condition (2. 2. 16b) we get 


w 


// (? + | Ui2 x^ r}:=r r w d©' dx* 


w 


w 


After a few simple manipulations (see Appendix A(i ) ) , 
we obtain 


I s = / [(A + B) 0(x ! , r ) - C | 0 x ,(x', r w )] dx' 

W - 00 p 

(2.3.7) 


where 



A = 

(r w /n)(r w - r) E(it/2, k)/[(a - b)a 1//2 ] 

(2.3.8a) 

B = 

[F( %/2, k) Ei.%/2, k.)] /(2ita 3 ‘ //2 ) 

(2.3.8b) 

C = 

r w F(tc/ 2, k.)/(7ia 1//2 ) 

(2.3.8c) 


a = (x - x?) 2 + (r + r w ) 2 , b = 4rr w , k = (b/a)^ 2 

(2.3. 8d) 

and F and E are complete elliptic integrals of the first 
and second kind respectively. 

The surface integral over the outer boundary can there- 
fore be expressed as 

I s = 8 / |(A + B) 0{x* / r w ) - C |0 xr (xS r w )] dx’ (2.3.9) 

T — CD ■ ' 
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where 


6 =0 for free air case 

= 1 for wall interference case. 


Using Eqs. (2.3.4) to (2.3.9) in Green's identity 
(2.3.1) / we get the following integral equation for the pertur- 
bation velocity potential 


0(x, r) ~ 


JC 1 s 8 (x 3 ) dx 1 

2 f 77 " ~~2 27172 

4vc p o [ (x - x * ) +r] / 


oo 

+ 6 / [(A + B) 0(x* / 

- 00 


r w } - 


c f 0 X ,(X', 


w 


)] 


dx' 


/// f g (x * / r 1 ) r* d©‘ dr‘ dx' (2.3.10) 

Q 


The integral equation for the axial velocity component 
is obtained by differentiating Eq. (2.3.10) with respect to x 
(see Appendix A(ii)). 


u(x, r) = -a- / — s'(x') dx* 

4itp o [(x - x') + r ] ' 

oo 

+ 6 / (A + B + D) u (x ' / r w ) dx* 

- 00 * 


+ SU 'P , (x, r; x* , r' , ©' ) gCx*, r* ) 


r' d©' dr* dx* 


U L + 6U w + U F 


(2.3.11) 


where 


D = PJ-X*.) x„r x 1 

P r w “ r 
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and subscripts L, w and P stand respectively for the linear/ 
wall and field contributions to the perturbation velocity. 

Equation (2.3c 11) is the non-linear integral equation 
for the axial component of the perturbation velocity whose 
solution must be obtained numerically. The numerical scheme 
for its solution is presented In the next section. 


2.4. NUMERICAL SOLUTION 

2.4.1. Discretization of the Integral Equation 

In order to discretize Eq. (2.3.11), x-r plane is 
divided into rectangular panels (Figure 2.2a). Making a cons- 
tant source approximation for each panel element, the field 
integral may be approximated by the sum 


u F (x, r) = 


S S 
k. 1 



(x/ r) g^ 


(2.4.1) 


where the indices k and 1 run along the x and r directions 
respectively . 

Similarly the wall integral in Eq. (2.3.11) may be 
replaced by the sum 


u w (x, r) = | a w (x, r) u^, 


1 = 1 


max 


(2.4.2) 


The aerodynamic co-efficients a ? and a^ are defined as 


kl 


v 6 


a-* 


ap (x, r) = * / dx* f r' dr* / 

- _Jk /- 

‘I 


2 % 


3¥ 


P ' ^-6 “- r,-h ' o’ 8*' 


as 1 


x ^.+6 

a*(x, r) = / (A + B + D) dx' 

Xk-6 


(2.4.3) 


(2.4.4) 
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where 26 (=A x) and 2h respectively are the width and height 
of each panel element, see Figure 2.2b. The indicated integ- 
rations in Bqs 0 (2.4.3) and (2.4.4) are performed after 
approximating the values of the elliptic integrals occurring 
in the integrands by their values at the centre of the element. 
The expressions for a p and a.^ are given in Appendix B. 

The first term in Eq. (2.3,11) represents the linear 
velocity distribution for subsonic axisymmetric flow past 
slender bodies. The integral occurring in this term must in 
general be numerically evaluated. A careful procedure is 
necessary for evaluating this integral for points close to the 
surface where the integrand has a singular behaviour. We have,- 


u L (x, r) 


K r s’U) _ } s u (x t ) n 

2 £"r, ”72 2 -.1/2 “ f 7~, rrt 2 -,1/2 dx J 

4 itp [(x - 1 ) + r J ' o [ (x - x‘) + r j 


K r s'(l) 

47ip 2 [(x - l) 2 + r 2 ] 1 ^ 2 


+ s" (x) In 


x - 1 + [(x - 1) 


, , 2. , l. 

x + (x + r ) 


z zn 

± r I 

2 , 1/2 




+ / r . dx 1 ] (2.4.5) 

o [(* r x ') 2 + r 2 ] 1/2 

where s*(o) has been assumed to be zero. The integral occurring 
in the third term of Eg. (2.4.5) is well behaved for points 
close to or on the surface and can be evaluated by quadrature. 


2.4*2. The Difference Operators 

In the discretization scheme to be adopted, we shall 
be using first and second order accurate difference operators 
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to evaluate the derivatives. They are defined as follows. 

Let the grid lines parallel to the y-axis be spaced a 
uniform distance A x apart in a Cartesian mesh and a function 
J be defined at the grid points/ see Figure 2.2c. A second 
order accurate central difference operator on such a mesh may 
be defined as 

C 3 o 25 J ]i = (J i + l - < 2 - 4 

where the subscript c stands for central differencing. 

Similarly other first order and second order accurate 
backward and forward difference operators may be defined as 


[3' 2) 

J ^i " 

(3J a - 4J i _ 1 + J ± )/2 Ax 

(2.4.7) 

[3f 2) 

J ^i " 

(- 3CT + 4J i+1 - J. +2 )/2 AX 

(2.4.8) 

[a cu 

^i = 

(J i ' J i-1 )/Ax 

(2.4.9) 

[ai 1J 

J h = 

" (J i ’ J ± + 1>/ Ax 

(2.4.10) 


where the subscripts b and f denote backward. and forward diff- 
erencing respectively, 

2.4.3. The Non-linear Source Term 

The function g(x, r) appearing in Eq. (2.2.17) represents 
the field source distribution. Different computational schemes 
may be formulated depending upon the form in which this term is 
written. In the present work three different representations 
for the evaluation g were tried. Denoting the schemes as Si/ S2 
and S3 we write/ 
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3 2 

g(x, r) « (j-) ” e 3 - [f(u)u x ] *. SI (2.4.11a) 

= uu x - e ^ [f (u)u x ] : S2 (2.4.11b) 

= uu^ - e f ^ u ) u xx : s3 (2.4.11c) 

In [73]/ Eq. (2.4.11b) has been called a quasi- 
conservative scheme and Eq. (2.4.11a)/ a fully conservative 
scheme. Schemes SI/ S2 and S3 as formulated above/ are comple- 
tely equivalent for purely subsonic flows because the artificial 
viscosity term vanishes identically, while differencing of the 
inviscid term in Eqs. (2.4.11) in conservative and non-conser- 
vative form are both equivalent in the sense that they converge 
to the same solution in the limit as the mesh width goes to 
zero. 

For super-critical flows with shock, the three schemes 
are capable of giving solutions which differ considerably in 
the shock region. This is so because the captured shock jumps 
are functions of the form in which the governing equations are 
differenced/ the correct jump being the one obtained when the 
governing equations are differenced in a form expressing a 
conservation law. Since the inviscid equation in the present 
formulation is already written in the quasi-linear form, schemes 
SI and S2, in which the artificial viscosity term is in conser- 
vative form, may be said to correspond to the quasi-conservative 
formulation while scheme S3, with the non-conservative viscosity 
term, may be said to correspond to the non-conservative quasi- 
linear formulation. 
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Introducing a switching function p in the artificial 
viscosity term 


H = 0 


u < 1 


= 1 u > 1 


(2.4.12) 


where u = 1 corresponds to sonic conditions , the field source 
term may be computed from 


2 

g X j = [3c 2) (f-)] ^[3b X) (M(u)VJij 5 S1 (2.4.13a) 

i j 

= [uu] - x (p f(u)u ] : S2 (2.4.13b) 

x ij D x 1J 

= [uu,] - X f ij C3b X) V ] ij ! S3 (2.4.13c) 

where X =e/Ax and 

u^j = [ u ^ij for J-uterior points 

(2.4.14) 

= [ 9^ ^ u]^. or [ 9 f 2 ^ u 3±j for boundary points 

The numerical viscosity parameter X helps in regulating the 
magnitude of the dissipation term and hence the stability of 
the iteration scheme. Use of Jameson's switching function p 
in conservative form in the artificial viscosity term leads to 
the correct sonic and shock point operators used first by 
Murman [6] in his conservative scheme. 

In the original non-conservative scheme of Murman and 
Cole [5], use of type dependent differencing for the x-deriva- 
tives led to the use of an elliptic operator at the first sub- 
sonic mesh point downstream of the shock (to be referred to 
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hereafter as the shock point). This was later called the non- 
conservative scheme by Murman [6]. Since non-conservative 
results are often closer to experimental results,, they are of 
some interest, although the shock jumps calculated by such 
schemes do not satisfy the conservation laws. Hence in the 
present work results have also been obtained by using the elli- 
ptic operator at the shock point. The field source term, at 
the shock point is accordingly calculated from 

2 

iij = C a4 2) (j->] ij s S1 (2.4.15a) 

= u. . [dl u], . i S2 and S3 (2.4.15b) 


2.4.4. Iterative Techniques 

The discretized approximation to the integral equation 
may be summarized as follows: 



u 


Li j 



(2.4.16) 


u_ . . = 2 2 a " . . g, , 

Fij ^ ^ Fxj y kl 

^Wij = J ^Vij u kl / 1 1 max 
k 


(2.4.17) 


(2.4.18) 


The field source term g may be calculated from anyone of Eqs. 
(2.4*13), depending upon the scheme chosen. At the shock point, 
g_ is obtained using either Jameson's shock operator (referred 
to as JA) from (2.4.13) or the elliptic operator (referred to 
as EL) .from Eq. (2.4.15). 
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Equation (2.4.16) represents a system of non-linear 
algebraic equations which must be solved iteratively. Numerous 
iterative schemes are available for the solution of such 
systems [94], In the present case, this system of equations is 
solved by a Direct Iteration Scheme (DIS) as proposed here or 
by a Quasi -Newton Scheme (QNS) as was done by Piers and Slooff 
[733. The latter scheme/ with a super-linear convergence rafce, 
is generally faster than the direct scheme but is found to be 
unstable for the wall interference case. The direct scheme, on 
the other hand, is found to be quite inefficient for free air 
calculations. 

For the free air case, either g or u may be treated as 
the iterative unknown. For the wall interference case however, 
only u is suitable as the iterative unknown because the wall 
integral requires the value of u for its evaluation. Further, 
numerical experiments show that it is desirable to add terms 
representing artificial time derivatives of u and u to the 

X 

iterative formulation of Eqs. (2.4.16) to improve the stability 
and convergence properties of the scheme. Assuming that an 
initial estimate g^ is available at the n-th iteration, Eq. 
(2.4.16) may be rewritten in the following form for the free 
air case. 


(n) 
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“lii + u 


(n) 

Fij 


(2.4.19) 


R 


(n) 

ij 


- 5*. - a,(u‘ n) - 

y ij l v ij ij ' 


t (n) 

CE 0 (UJ. 

2 x ±j 


x ij 


(2.4.20) 



51 


Here g* is obtained from. Eq. (2.4.13) using the estimates u^ n ^ 

/ v 

and u^. 1 . and a 2 are iteration parameters used to control 

the magnitude of the time damping terms and hence the stability 
and convergence rate of the iterative scheme. The next iterate 
is obtained from 


5 (n. +1 ) 



OA 


:(rt) 


(2.4.21) 


where to is a damping factor for the correction vector. The 
value of to is usually less than unity and is different at sub- 
sonic and supersonic points. 

The correction vector Ag^ n ^ is obtained by solving a 
system of linear equations in the quasi-Newton scheme. In the 
direct iteration scheme, it is simply the negative of the 
residual. Thus 


Quasi -Newton Scheme: 


3tl(n) -(n) 
c ij A g kl 


- R (n) 

R ij 


Direct Iteration Scheme: Ag^ 1 ^ = 

AJ 


R (n) 

R ij 


(2.4.22) 

(2.4.23) 


The matrix C is computed from 


c kl(n) 




(2.4.24) 


The matrix c is approximated by its tri -diagonal elements and 
the correction vector is easily determined by an explicit inver- 
sion formula. 


For the wall case, assuming an initial estimate u 


(n) 


we have. 


i Urn 


■ am 







52 



The correction vector is calculated by a procedure similar to 
that outlined for Ag in the free air case. 


2.5, RESULTS AND DISCUSSION 

The numerical schemes were coded in FORTRAN and run on 
the DEC 1090 computer . Most of the calculations were carried 
out interactively on the time-sharing mode. In the following 
sections, results obtained for axisymmetric flow past slender, 
general parabolic bodies of resolution in free air as well as 
in the presence of perforated wind tunnel walls will be 
presented. 

The equation for the general parabolic body of revo- 
lution is given in Appendix C. The linear solution u^, given 
by Eq. (2.4.5), was computed numerically, evaluating the 
integral by Simpson' s rule . For the case of a body with 
maximum thickness at the centre, the numerically computed 
values were compared with exact values obtained by evaluating 
the integral in Eq. (2,4.5) analytically. The two results were 
found to be in agreement within an error margin of about 0.5 
percent . The linear solution was then used as an input for the 
iterative, solution, of the non-linear problem. 
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Finally it is remarked that the' value of the artificial 
viscosity parameter X was set equal to unity in all the foll- 
owing calculations except when stated otherwise. 

2.5.1. Grid Effects 

Following normal practice, the grid related parameters 
for obtaining satisf actory solutions were established by numer- 
ical experiments. These results are discussed first. 

Figures 2.3 (a to e) present the effects of the grid 
related variables on the surface pressure distribution for a 
parabolic body of revolution without sting and with maximum 
thickness section located at the mid-point of the axis. The 
free-stream Mach number and body fineness ratio chosen to study 
these effects were 0.98 and 10 respectively. The numerical 
calculations were carried out for a typical free air case" using 
scheme and shock operator JA. 

Figures 2,3a and b show respectively the effect of 
axial and lateral grid extents. For the three grids used in 
each of these figures (namely j x i * 10 x 40, 10 x 35 and 
10 x 31 in Figure 2,3a and 10 x 35, 7 x 35 and 6 x 35 in Figure 
2.3b), it is seen that the surface pressure distribution obtained 
on the smallest grid differs hardly from that obtained oh the 
largest. Hence it is concluded that a computational domain 
bounded by -0.20 < x < 1.20 and 0 < r < 0.35 i« sufficient 
to obtain satisfactory results on the body surface. This may be 
compared with the axial grid extent of -2.36 < x < 3.36 used by 

•alley [84j in his finite-difference calculations. For = 













59 


0.99, the above range of r on the transformed plane corresponds 
to 0 < r _< 2.5 on the physical plane • 

Further economy of grid points was achieved by progre- 
ssively reducing the number of streamwise grid points at every 
j-level, moving away from the body axis (see Figure 2.3c). The 
figure shows that with I = 2, the number of grid points could 
be reduced to 161 from 245, leaving the computed pressure 
distribution on the surface unaffected. The CPU time gets 
nearly halved with their reduction in the number of grid points. 
Calculations presented here were however carried out with 1=1 

A 

unless stated otherwise. Further it is pointed out that for 
the wall interference calculations all the grid points on the 
wall boundary are esential for the evaluation of the wall 
integral and hence I = 0 for all these calculations. 

A 

Figure 2.3d shows the effect of mesh size on the surface 
pressure distribution. As expected, a coarse mesh causes a 
slight smearing of the shock. This is evident from Figure 2.3d 
for Ax = 0.05. Figure 2.3e shows the effect of a stretching 
factor in the r-direction on the surface pressure distribution. 
This factor increases the aspect ratio of the rectangular cells 
in the r-direction as one moves away from the body axis* Thus 
for a fixed lateral grid extent fewer number of j-levels are 
required than if a square mesh were employed. For the two 
values of SF = 1.2 and 1.4 vised in Figure 2.3e, the computed 
pressure distributions turn out to be identical. It might seem 
that still larger values for the stretching factor could be 
used to reduce the number of grid points further, but since the 
iteration parameters go , a 1 and a 2 are related to mesh aspect 
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ratio,- a larger SF could, and often did, lead to divergence of 
the iteration scheme. For example with M = 0.99 and' PR = 10, 

CD 

it was not possible to obtain convergence for the free air 
case with SF > 1.0, whatever the values of o , <x^ and used. 

The values of SF used in the present calculations generally 
varied from 1 to 1.2 depending upon the free stream Mach number 
and body fineness ratio, the smaller value being used for the 
stronger super-critical cases. 

2.5.2. Comparison of the Results of the Schemas SI, S2 and 

S3 and Shock Operators JA and EL 

Figures 2.4 and 2.5 present a comparison of the computed 
surface pressure distributions for a parabolic body of resolu- 
tion without sting in free air, using the schemes SI, S2 and S3. 
Results obtained using both the Jameson (JA) and the elliptic 
(EL) shock operator are presented in the figures. Figure 2.4 
displays the pressure plots for =0.99, FR = 10 whereas 
Figure 2.5 shows the same plots for = 0.98, FR = 6, 

Considering Figure 2*4a, we see that scheme S2 shows 
an over expansion just ahead of the shock while elsewhere, the 
three results are in close agreement with one another. It is 
noted however that scheme SI needs just about half as many 
iterations to converge as do the other two schemes. The shock 
operator used in this case was JA. Figure 2.4b shows the 
results of the same calculation using the EL shock operator. 

The scheme S3 in this case went into an oscillation near conver- 
gence. Schemes SI and S2 converged in about 50 iterations. 

Here again the pressure distribution obtained from scheme S2 
is seen to overshoot just before the shock. 




Fig. 2.4. Comparison of the Schemes SI , S2 and S3 : Surface Pressure 
Distribution on a Parabolic Body of Revolution in Free Air . 
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Comparison of the Schemes SI, 52 and S3 : Surface Pressure 
Distribution on a Parabolic Body of Revolution in Free Air . 
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Figure 2.5 shows the results of the above calculations 
repeated with M = 0,98 and FR = 6. Comments made in relation 

CD 

to Figure 2,4 apply almost identically to the results presented 
in Figure 2,5. In addition it is noted that in Figure 2.5b, 
which shows the pressure distribution obtained with an elliptic 
operator, scheme SI also results in a slight over-expansion 
just ahead of the shock. The significant difference in the 
shock location calculated by scheme SI and that of S2 and S3 
for this case is due to the occurrence of this over-expansion. 
It was thought that such an over-expansion occurred because of 
insufficient dissipation. The calculations were therefore 
repeated for higher values of X but it was found that the 
non- uniformity persisted and the shock location remained unch- 
anged. Here again it is noted that scheme S2 in Figure 2.5a 
and S3 in Figure 2.5b went into oscillations near convergence. 

In all these cases, the surface pressure distribution 
shows a smooth development over the subsonic region ahead of 
the shock, including the region around the sonic point. The 
subsonic portion is, as expected, insensitive to the scheme 
used. There is also a small region of smooth deceleration, 
instead of an abrupt transition from the point of the velocity 
peak to the beginning of what can be regarded as the shock 
transition region. 

Figures 2.6a and b show comparisons of the surface 
pressure distributions obtained on the same body in free air, 
using the elliptic and Jameson shock operators, for = 0.99, 
FR « 10 and M =0.98 and FR = 6 respectively. The SI scheme 

CD 

was used for these calculations. In both cases it is seen that 




Fig. 2.6. Comparison of the Shock Operators (SO) : Surface Pressure 
Distribution over a Parabolic Body of Revolution . SCH:S1» 
Free Air . grid : 7x 35 . 
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the elliptic shock operator produces a weaker shock, located 
about half a mesh width ahead of that predicted by the Jameson 
operator. For weaker super-critical flows however, this effect 
does not emerge clearly since in such cases, the present method 
tends to smear out the captured shocks. 

From the above discussion it is clear that scheme SI, 
together with the Jameson or elliptic shock operator is the 
best choice from the point of view of computing time, as well 
as, obtaining satisfactory solutions. Results of further calcu- 
lations presented here were therefore obtained using the SI 
scheme together with either the Jameson (JA) or elliptic (EL) 
shock operator. 

2.5.3. Comparison of Present Results with Other Experimental 
and Numerical Results 

A comparison of the results of the present computations 
with other theoretical and experimental results is presented in 
this section. 

(a) Free Air Results 

Figure 2.7 shows the computed results in free air for 
a parabolic body of revolution with sting at x = 0,854, FR = 10 
and at a free-stream Mach number 0,975, using the EL operator. 
Pressure distributions on the body surface and at four radial 
locations away from the body are compared with the experimental 
results of McDevitt and Taylor [95], The agreement is seen to 
be reasonably good. Figure 2,8 shows the pressure distributions 
obtained on general parabolic bodies of revolution with axial 
locations of maximum body thickness at 0*3, 0.4, 0.6 and 0.7 




» t ' 

27. Comparison of Pressure Distribution for Flow 
over a Parabolic Body of Revolution at Different 
Radial Locations with Experimental Data. 

M w s 0.975 , FR s 10 , SCH : SI ( SO : EL , Free Air . 
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and PR «= 12. The free- stream Mach number is 0.975. The agree- 
ment with experimental data [96] is again in general satisfac- 
tory. The discrepancy near the body sting may be caused by 
the failure of small disturbance theory in this region. Use 
of Jameson shock operator gave almost the same results for 
these cases since the shocks are comparatively weak and smeared 
for all of them. 

Figures 2.9a and b show comparison of the present 
results with results of finite-difference calculations. Present 
results, obtained using the JA shock operator, are plotted for 
a body with sting in free air. In Figure 2.9a with M = 0.99 
and FR = 10, it is seen that the non-conservative shock loca- 
tion from Bailey 1 s [84] finite-difference calculations, is about 
half a mesh width ahead of the location predicted by the present 
quasi-conservative calculations. The pressure distribution 
obtained by using the EL operator for this case is plotted in 
the same figure and is seen to be closer to Bailey's results. 
Figure 2.9b shows the results of the same calculation for 
M = 0.98 and FR = 6. This is compared with Krupp's [97], 
non-conservative finite-difference calculations. As in the 
previous case, here too the shock location predicted by Krupp* s 
scheme is about half a mesh width ahead of the present shock. 

In both figures it is noted that the surface pressure distribu- 
tion calculated by the present scheme, in the fore region of 
the body is somewhat flatter compared to the finite-difference 
results. But for this discrepancy, the results of the present 
calculations are in fairly good agreement with experiment as 
well as with finite-difference results. 
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Fig. 2.9. Comparison of Computed Pressure Distribution for Parabolic 
Body of Revolution with Finite - Difference Results , Free Air . 
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( b ) Wall Interference Resul ts 

Figures 2.10 and 2.11 show, from calculations using an 
EL operator , the effect of the presence of a perforated wall 
on the surface pressure distribution for bodies of two different 
fineness ratios at two different free stream Mach numbers. 

In general it is seen that the presence of the wall 
has a significant effect on the computed surface pressure 
distributions. For the weaker super-critical cases, with = 
0.975, PR = 10 or M = 0.94, FR = 6 (Figures 2.10a and 2.11a), 
the effect is more pronounced in the shock region and relatively 
small over the subsonic and accelerating supersonic regions of 
the flow. For the stronger super- critical cases with M^ = 0.99, 
FR = 10 and M = 0.98 and FR = 6, the surface pressure distri- 
bution changes over a greater region including the accelerating 
supersonic portion and the region of shock transition. Figures 
2.10b and 2.11b. It is seen that increasing porosity causes 
the shock to move upstream by about 2 percent of the chord. For 
smaller free- stream Mach numbers, with fixed FR, there is a 
tendency for the shock to get smeared with increasing porosity. 

In general the surface velocities attained are smaller as the 
porosity factor increases. This makes physical sense because 
if part of the mass flow is bled away, the velocities must be 
smaller. The solid wall case with p = 0 gives the maximum 
surface velocities. However convergence could not be achieved 
for this case for the stronger super-critical cases (for example, 
= 0.99, FR = 10) . 



Effect of Porous Wall Boundary on the Surface 
Pressure Distribution for a Parabolic Body of 
Revolution with Sting. 
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over a Parabolic Body of Revolution . SCH : 51 , SO: EL. 
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Figures 2.10 and 2.11 also show the experimental data 
obtained for these cases plotted along with the present theore- 
tical results. The comparisons shown by Figures 2.10a and b 
suggest that the experimental shock may be approached by 
increasing porosity. 

Figures 2.11a and b show the results of the present 
calculations compared with experimental data [98] for a body 
with FR = 6 . Here the agreement with free air results seems to 
be better than that with wall interference. As for the shock, 
the experimental shock strength seems to be in good agreement 
with the calculated value for M =0.94 whereas for M =0.98, 

OD QD 

the calculated shock is somewhat stronger. The shock locations 
in both cases are however in agreement with experiment. 

2.5.4. Convergence Characteristics of the Iteration Scheme 

As mentioned in Section 2.4.4, free air calculations 
were done by iterating on g and wall interference calculations 
by iterating on u. The Quasi- Newton Scheme (QNS) or the Direct 
Iteration Scheme (DIS) could be employed for both. QNS proved 
to be more efficient for the free air case whereas DIS was 
found superior for the wall interference calculations. 

Figure 2.12 shows, for a fixed w, the effect of the 
iteration parameters and on the convergence history for 
the QNS . These calculations were done in free air for a body 
with no sting, with FR = 10 and M ^ = 0.99. The value for o> 
is given in pairs, the first value being used at subsonic points 
and the second at supersonic points. Since we had no apriori 
estimates for oc^ and a^, the solution for a prescribed free 
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stream Mach number and thickness was attempted by setting 
a-j_ = <*2 = 0- w was set at some value usually less than unity , 
the value used at supersonic points being 0.1 less than that at 
sufo«*s©imie points. A larger value for e generally implies 
faster convergence but could also lead to divergence depending 
upon the free stream Mach number and thickness. 

Figure 2.12 shows , for the QNS, | j R | i of the residual 
plotted against the iteration number. The two curves demons- 
trate the effectiveness of the artificial time-dependent terms 
in achieving faster convergence. The number of iterations 
needed with ce^ =7.0 and = 0.05 is 48 whereas with = ct^ = 0, 
76; iteration were needed for the same. For this case liRi i 
showed a more or less uniform decreasing behaviour from the 
maximum value at the beginning of the iteration to its value at 
convergence. The oscillatory behaviour noticed for the plane, 
two-dimensional case (see Section 2.5.5) and in the wall case 
below, was not evident in the convergence performance of QNS 
for the free air calculations. 

Figures 2.13a and b show the development of the corres- 
ponding surface pressure distributions for the two cases discu- 
ssed above. It is seen from Figure 2*13a that with = a 2 ~ 
the subsonic points attain their converged value (to plotting 
accuracy) in as few as 10 iterations. The rest of the iterations 
serve to adjust the location and strength of the shock. Even 
so, the surface pressure distribution practically stabilises in 
about 30 iterations. The convergence criterion for the present 
calculation being rather stiff (about 0.1 percent on ) , the 








Fig. 2.13. iterative Development : Surface Pressure Distribution on a 


Parabolic Body of Revolution in Free Air. 
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iterations were carried out for longer than was necessary to 
realise plotting accuracy which was about 1 percent. 

Figure 2.13b shows the iterative development of the 
solution with ck-^ =7.0 and ~ 0.05. Introduction of the 
artificial time-dependent terms is seen to result in a practi- 
cally converged solution being achieved in about 10 iterations. 
These figures also show that no large amplitude fluctuations 
develop in the solutions even for this fairly strong super- 
critical case,, neither in the subsonic region nor in the super- 
sonic or in the shock regions. This indicates the possibility 
of attempting non-linear extrapolation of the solution to the 
converged solution after a sufficient number of iterations in 
the initial stages have been carried out. 

Figure 2.14 shows similar results on the convergence 
characteristics of the Dis and QNS on u for the wall interference 
case. The free-stream Mach number is 0.99 and FR = 10. The 
porosity was set at p = 0.75. Four of the lines in Figure 2.14 
show the error norm plotted against the number of iterations 
for DIS. The fifth line traces the convergence history of QNS. 
The damping factor^ for this case was (0,2, 0.3). It is to be 
noted here that the second value in the paranthesis is larger 
than the first. This was the only case for which a smaller 
value of w for th- supersonic points and larger value for the 
subsonic points failed to give convergence with or without the 
art if icial tim'"'— do— indent terms. Even values smaller than (0.2, 
0.1) led to divergence of both DIS and QNS, whereas with 03 
(0.2, 0.3), convergence could be achieved albeit only with the 
presence of the artificial time-dependent terms. Figure 2.14 
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shows that, with ct^ = 1.5/ =0.10/ convergence could be 

achieved in just 35 iterations. Comparison of this case with 
the one for ~ 0 shows that convergence here is smooth 

and fast. The oscillations noticed in the behaviour of !| Rj| 
plot for a i “ = 0 are damped out rapidly with the addition 

of the time-dependent terms. This is also confirmed by a 
study of the surface pressure plots shown in Figures 2.15a 
and b. 

Figure 2.15a shows the development of the surface 
pressure distribution with a ]_ = a 2 ~ The a PP earance °f a 
prominent non- uniformity in the surface pressure may be noticed 
just ahead of the shock at the end of about 20 iterations. At 
the end of 40 iterations/ this non- uniformity has grown larger. 
The surface pressure thereafter tends to smoothen out as the 
iteration proceeds and the error norm correspondingly shows a 
decreasing trend. This cas.i was found to converge extremely 
slowly as the iterations were carried on beyond 60. Figure 
2.15b shows the development of the solution beginning with the 
initial solution. No non- 'uniformities appear in the solution 
as the iteration proceeds and convergence is apparently quite 
smooth. The solution at the end of 20 iterations is seen to 
be practically indistinguishable from the final converged 
solution which was achieved in 35 iterations. 

The QNS also produced converged solutions for this case 
in about 90 iterations. The convergence history for this case 
is also displayed in Figure 2.14 and is seen to be inferior to 
that of the direct iteration scheme . This conclusion was also 
borne out by tests carried out for different values of the 
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Fig. 2.18. Iterative Development: Surface Pressure Distribution on a 
Parabolic Airfoil. 0.85 r = 0.10 grid: 10x30 
(0.4, 0.3) a, = 0.0 a 2 = 0.0 
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Fig. 2.19, Iterative Development : Surface Pressure 
Distribution on a Parabolic Airfoil . 
M«>=G.85 t - 0.10 grid : 10x30 
w = (0.5, 0.4) a, = 3.1 a 2 = 0.2 



85 


flow parameters as well as various combinations of the values 
of w , and 

Admittedly this study of the artificial time-dependent 
damping terms merely demonst.pat.es their utility in imparting 
stability to the iterative scheme. We have found yet no way 
of arriving at apriori estimates for the parameters to t and 
a 2 » The usefulness of these terms, in the absence of such 
estimates is restricted to achieving convergence for the stronger 
super-critical cases which might otherwise diverge. By a trial 
and error procedure, these cases may be made to converge by 
choosing appropriate values for to, and a 2 * 

2.5.5. Some Results for Plane, Two-dimensional Symmetric 

Flows 

A code analogous to the axisymmetric code was also 
developed for plane symmetric flow past airfoils. The results 
of Piers and Slooff [73] were in the main reproduced with this 
code. The results obtained using scheme SI together with JA 
operator was found to give the best agreement with those of 
finite-difference calculations, as shown in Figure 2.16 for a 
10 percent thick parabolic arc airfoil at = 0.825 and 0.85. 
Numerical experiments were also carried out for this case by 
including artificial time-dependent terms. Reduction in the 
iteration counts comparable to those reported above for the 
axisymmetric case, were observed for plane flows also. These 
observations are summarized in Figures 2.17 to 2.19. 



CHAPTER III 


PLANE, TWO-DIMENSIONAL TRANSONIC FLOW PAST 
AIRFOILS - FULL-POTENTIAL SOLUTIONS 

3.1. INTRODUCTION 

In this chapter, we present a hybrid calculation 
procedure, combining the Integral Equation Method with elements 
of the finite-difference technique to compute full-potential 
solutions to two-dimensional transonic flows over airfoils with 
large embedded supersonic regions and weak shock discontinuities. 
In this approach, two simultaneous integral equations for the 
velocity components are formulated by considering the non-linear 
terms in the full -potential equation as a field source distri- 
bution, together with a distribution of sources and vortices on 
the mean line of the airfoil. 

To be able to compute super-critical flows with shocks, 
an artificial viscosity term is added to the field source term 
in the supersonic regions so that expansion shocks are excluded 
from the computed solutions while the compression shocks are 
replaced by continuous regions with steep velocity gradients. 

Since the integral equations for the velocity conponents are 
coupled and non-linear, an iterative procedure is necessary for 
the numerical solution of these equations. The non-linearity 
consists in the presence of the field source term which is a 
function of the velocities and their derivatives. In order to 
calculate these derivatives, a body conforming mesh is generated 
through a sequence of co-ordinate transformations. The derivative 
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calculations are carried out on the transformed plane and the 
derivatives on the physical plane are obtained using the 
Jacobian of the transformation. 

The iterative procedure starts by assuming the field 
source distribution to identically vanish and computes the 
incompressible solution to satisfy the tangency boundary condi- 
tion on the body surface. The incompressible solution gives 
rise to a field source distribution which results in additional 
induced velocities/ disturbing the boundary condition on the 
body surface. The internal singularity distribution is accor- 
dingly redetermined in order to satisfy the surface boundary 
condition iteratively. The non-linear problem thus reduces at 
every iteration to the solution of a Poisson problem. The 
procedure is repeated till convergence is achieved in the 
entire flow field. 

The chapter concludes with a discussion of the computed 
solutions for sub-critical and super-critical flows past a few 
selected airfoils. 

3.2. PROBLEM FORMULATION 
3.2.1. The Inviscid Equations 

Consider plane, two-dimensional flow past an arbitrary 

airfoil whose chord is aligned with the x~axis and the leading 

edge is located at the origin of a Cartesian c o -ordinate system. 
Figure 3.1. 

The equations governing compressible, inviscid and 
irrotational flow of a perfect gas can be written as 



xCosa ♦ y Sin Ct 



Mathematical Model of Plane, Two-Dimensional, Compressible Flow 
Past an Airfoil. 
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Mass conservation ; ( PVg ) = 

„ ^ . 2 Y - 1 rs2 

Energy conservation; o + — - Q 

where q is the velocity potential, 
speed of sound and Q the fluid speed* 

For isentropic flows in which energy is conserved, it 
is not possible to conserve both mass and momentum across a 
shock. For an approximation in which mass is sought to be 
conserved, the momentum deficiency then provides an estimate 
of the wave drag. 

To the conservation relations (3.2.1) and (3.2.2), the 
following isentropic relations may be added; 

y 22 Y— 1 

P = P and c = P* (3.2.3) 

In Eqs. (3.2.1) to (3.2.3), all velocities are made 
dimensionless by the free-stream velocity while pressure and 
density are made dimensionless by their respective free~stream 
values. The length scales are measured with the airfoil chord 
taken as unity. 

Combining Eqs. (3.2.1) to (3.2.3), the so called gas 
dynamic equation is obtained 

fc2 ' u2) ^xx - 2UY 5xy + (o 2 - V 2 ) § yy = 0 (3.2.4) 

where U and V are the velocity components in the x and y 

2 2 2 
directions and Q = U + V . 

Equation (3.2.4) may also be written in the form 

5 + § ss [u JL (^—) + v — (~)] s a (3 2 5) 

* xx j yy 2 LU 9x v 2 ; 3y v 2 }1 g 

c 


0 


(3.2.1) 


1 . Y - 1 

“T-+ — 

ao 


(3.2.2) 

P the density, c the local 
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where g, as a function of the velocities and their derivatives, 
represents the non-linear terms on the right acting as a field 
source distribution. 

The total potential § may be split into a free-stream 
component and a perturbation component. Thus, 

5 (x, y) = x cos a + y sin a + 0(x, y) (3.2.6) 

where 0 is the perturbation potential and a is the angle of 
attack. Accordingly the velocities are given by 

U = cos a + u, V = sin a + v (3.2.7) 


where u and v denote the perturbation velocity components in 
the x and y directions respectively. 

In terms of 0, Eq. (3.2.5) simply becomes 

^XX + ^yy “9 (3.2.8) 

Let the airfoil shape be prescribed as 


Y 


U/l 



0 ^ x b l 


where the subscripts u and 1 stand for the upper and lower 
surfaces respectively. For a solid surface, the body boundary 
condition is specified by setting the normal component of the 
flow velocity to zero on the body surface. In terms of the 
perturbation velocity components this can be expressed as 


p u,l (x > u 


cos a F’ . (x) - sin a 

U./ JL 


on 


U,1 


Vl w 


(3.2.9) 
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where prims denotes differentiation with respect to x. 

Since the perturbation velocity components must vanish 
at infinity we have 


_ , 2 2 x l/2 

u, v -» 0 as (x + y ) / .» od 

In addition to Sqs , (3,2,9) and (3.2.10) , the 
condition of smooth flow at the trailing edge must be 
The local pressure co-efficient is calculated 


(3.2.10) 

Kutta 

satisfied. 

from 


P 







(3.2.11) 


3.2.2, Super- critical Hows, Artificial Viscosity 

For super-critical flows, as in the TSPE formulations, 
it is necessary to add an artificial viscosity to the governing 
equation in order to exclude expansion shocks and capture only 
the compression shocks. To arrive at an appropriate artificial 
viscosity model for the integral equation method, it will be 
instructive to study the way in which artificial viscosity is 
introduced in finite-difference calculation schemes [25]. 

The basic idea behind Murman’s [ 5 ] mixed differencing 
scheme is to introduce directional bias as well as to incorpo- 
rate the correct domain of dependence in the difference scheme 
at supersonic points. In the small disturbance formulation, 
where the flow direction is assumed to be more or less aligned 
with the x-co-crdinato axis, this is accomplished by simple 
upwind differencing of the x-derivatives at the supersonic 
points (see Section 2 ,2 .4 ) „ For the full-potential formulation 
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however , the upwind direction is not known in advance. ever- 

theless, for flows with supersonic zones of moderate size/ 

satisfactory results have been obtained by the use of a simple 

upwind differencing scheme in which central difference formulae 

are used at subsonic points for all derivatives e nc ^ upwind 

formulas for 5 and § at the supersonic points. The upwind 
xx xy 

differencing of Eq. (3„2.4) at the supersonic points then intro- 
duces an effective artificial viscosity 

Ax [(U 2 - c 2 ) U + UW ] 

When the flow is not perfectly aligned with the x-axis , 

2 2 2 2 one 

there exist supersonic points at which U<c< u_f * ne 

characteristic lies ahead of the y-co-ordinate li ne at suctl a 

point, so that the difference scheme does not have the correct 

region of dependence, see Figure 3.2a. Correspondingly, the 

artificial viscosity Ax(U 2 - c 2 )U , introduced by the use of 

the upwind formula at such a point becomes negative. Despite 

this fact, schemes of this type have proved satisfactory for 

flows with moderate supersonic zones. 

The treatment of flows with large supersonic zones in a 
curvilinear co-ordinate system, suited to the geometry of the 
problem, requires the use of a more elaborate dif f erence scheme 
in which the direction of upwind differencing i s independent of 
the co-ordinate system and is rotated to conform to the local 
flow direction. Setting up a local Cartesian co-ordinate system 
and denoting the local streamwise and normal directions by the 
s and n co-ordinates, the quasi-linear equation (3.2.4) can be 


written as 
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Fig. 3.2a. Simple Upwind Difference Scheme j [25] 



0 


(3.2.12) 


2 2 

Cc - cr ) 


ss 


, 2 
+ c 5 


nn 


Since U/Q and V/Q are the local direction cosines, § 

ss 

and § nn may be expressed in terms of the original co-ordinates 
as 


ij (u 2 *^ + 2UV 5xy + v 2 Syy ) (3.2.13) 

^2 < v2 Sxk * 2UV 5 xy +u2 5yy> t 3 - 2 - 14 ) 

Then at subsonic points, central difference formulas are used 
for both 3 ss and 5 nn » At supersonic points, § is calculated 
using central differences, but upwind formulas are used for 
evaluating all the derivatives occurring on the right of Eq. 
(3.2.13). Such a scheme consequently has the correct region of 
dependence at the supersonic points as shown in Figure 3.2b. 

Upwind differencing in this case introduces an effective artifi- 
cial viscosity 

2 

(1 - -2-5 [ Ax(U 2 U + UVV ) + Ay(UVU + V 2 V )] (3.2.15) 

q2 u xx xx - 1 yy yy . v 

which is symmetric in x and y. 

Ip the integral equation formulations, the artificial 
viscosity term is to be added explicitly to the field source 
term. This objective is realised by modelling the added arti- 
ficial viscosity term to correspond to that given by Eq. (3.2.15). 
Going back to Eq, (3,2,12), we see that if at the supersonic 
points 5 is represented by upwind difference formulas, this 
is equivalent to the addition of an artificial viscosity 


ss 


§ 


nn 
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= - As 


2 

P Q 5 


sss 


(3.2.16) 


where p is a switching function defined as 

r c-2-, 

pi = max [ 0 , 1 - -~ 5 * ] 

q 


(3.2.17) 


Since the mass continuity equation is equivalent to 
multiplying Eq. (3.2.4) by P /c? g the conservative form for the 
artificial viscosity term may be approximated as 



(3.2.18) 


It may be verified that the highest derivatives in the 
leading term of Eq. (3.2.18) are the same as those in Eq. 
(3.2.15), see Appendix D. 

Equation (3.2.18) is written in terms of the streamwise 
derivatives and hence can be applied conveniently on the surface 
grid points. For the interior points however, it is necessary 
to write the artificial viscosity in terms of the x and y deri- 
vatives as (see Appendix P) . 

T = 9x + "dy (3.2.19) 


where the artificial viscosity functions G and H are defined as 

,2 


Ax t- ** u h <f-> 

c 


(3.2.20) 


H « 


Ay -o v v 


JL /SC 

3 y ^2 


(3.2.21) 


With the addition of the artificial viscosity term to the right 
of Eq. (3.2.5), the modified field source term may be written as 
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g = g (x, y) + T (3.2.22) 

where T is given by either (3.2.18) or (3.2.19). 

In Eq. (3.2.22) it is noted that while the inviscid 
term g is not in conservation form, the artificial viscosity 
term is in conservation form. Hence a scheme based on Eq. 
(3.2.22) may be said to correspond to the quasi -conservative 
scheme of Bauer et al [ 99 ] . 

3.2.3. Artificial Compressibility 

The fact that the inviscid field source term may also 

3 P 

be expressed as g = - Q/p 5 — , suggests that the idea of arti- 
ficial compressibility may be just as easily implemented in the 
present scheme as artificial viscosity. The relative merits of 
the two approaches in a numerical procedure are not clear yet 
although Harten[lOO] has claimed that the artificial compressi- 
bility method is superior to artificial viscosity, in particular 
for capt Taring contact discontinuities. 

A simple approach to introducing artificial compressi- 
bility in the inviscid equations is to modify the continuity 
equation as 

3 3 

ax 5 x^ + (TT 5 ) = 0 (3.2.23) 

where the modified density is given by 

P = p - n As P s (3.2.24) 

Equation (3.2.24) may be arrived at as follows. 

The artificial viscosity function given by Eq. (3.2.20) 
is 



98 


G 


P 

- Ax p. U 


3 

3x 



. 3P 


since 


3 P 

$x 


Similarly 


H 


Ay 




v 


3 P 

ay 


_ 3 ^ 
c 2 3x 



) 


The form of these two functions immediately suggests 
writing the artificial compressibility term in the form given 
by Eq. (3.2.24), 

The modified field source term g for the artificial 
compressibility method is then 


q3P , 

g “ " jr 3 S ^ 3 * 2 “ 25 ) 

with P defined by Eq. (3.2.24). 

Other artificial compressibility terms may be found 
discussed in Ref. [lOl] . 


3.3. THE INTEGRAL EQUATION 

For the purpose of formulating the integral equation 
for the perturbation potential/ Eq. (3.2.5) may be regarded as 
a Poisson equation with a field source term g given by Eqs. 
(3.2.22) or (3.2.25). The solution for 0 may then be expressed 
as the sum of a general solution satisfying the homogeneous part 
of the equation and a particular solution. 

Since the general solution satisfies the Laplace equa- 
tion/ it may be built up by a suitable distribution of singu- 
larities. For lifting flows over two-dimensional airfoils / a 
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combination o£ sources and vortices are distributed on the mean 
line of the airfoil, as shown in Figure 3.3. The non-linear 
term g. can be thought of as a field source distribution. Hence 
the potential induced at a field point (x, y) is the sum of the 
potentials induced by the internal singularities and the; field 
sources. Assuming the source density to be a(x), vortex density 
Y(x) , we may write 


0(x, y) 


/ K s (x - x* , y - y (x* ) ) d(x* , y (x‘ ) ) dl 
+ / KyCx - x ? , y - y m (x l ) ) y (x* , y m (x* ) ) dl 


+ ff Kg (x - x' , y - y l ) g( x *, y* ) dx f dy* 

where the kernels are given by 


(3.3.1) 


Kg(x - x* / y - y* ) = In [ (x - x* ) 2 + (y - y* j 2 ] 1 / 2 

(3.3.2) 

Kv(x ” X’ , y - y l ) = tan” 1 (3.3.3) 

In Eq . (3.3.1), D is the entire plane excluding the area 

enclosed by the airfoil contour. The line integrals are over 

the curve C which denotes the mean line y (x) on which the 

m 

singularities are distributed. A small gap 6 is maintained 
between the leading edge and the tip of the mean line so as to 
render the velocities finite at the leading edge. In Eq. (3.3.1) 
hence, the first two terms represent the potential induced by 
the source and vorticity distributions, while the last term 
represents the effects of compressibility. 



#(x,y) Field point 
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The perturbation potential given by Eq. (3.3.1) 
satisfies the differential equation (3.2.5) in the region D 
with the field source term replaced by g, as can be verified 
by direct differentiation. Hence 0 is a solution of Eq. (3.2.5) 
in D. 

The integral equation for the perturbation velocity 
components may be obtained by differentiating the expression 
for 0. For brevity we adopt the following functional notation 
and write 


0(x, y) 

- V (K S' V 


u(x, y) 

- V CK Sx- K VSc > 

(3.3.4) 

v(x, y) 

= 0-f < K Sy' *V 

(3.3.5) 


where the subscripts jc and y denote the respective partial 
derivatives of the kernels. 

Equations (3.3.4) and (3.3.5) constitute a system of 
two integral equations for the perturbation velocity components 
u and v, which must be solved by satisfying the tangency boun- 
dary condition Eq. (3.2.9) on the surface of the airfoil- The 
far -field boundary condition is implicitly satisfied by the 
integral equations as can be verified for the first two terms 
by inspection and, for the field integral by making appropriate 
assumptions about the asymptotic behaviour of g. 

The solution of the Eqs. (3.3.4) and (3.3.5) is to be 
obtained by numerical methods. 
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3.4. NUMERICAL SOLUTION 

A glance at Eqs. (3,3.1), (3,3.4) and (3.3.5) suggests 
that it may be advantageous to solve Eq„ (3.3.1) for 0 because 
it involves the computation of fewer number of influence 
co-efficients . However in a hybrid calculation procedure of 
the kind attempted here, accurate evaluation of velocities 
from the perturbation potential using finite-*difference deri- 
vatives would demand a rather close spacing of the grid lines 
adjacent to the airfoil as in finite-difference methods. On 
the other hand, dealing with Eqs. (3.3.4) and (3.3.5) directly, 
involving the computation of the influence co-efficients for 
both the velocity components, permits accurate evaluation of 
the same, especially at the surface points. This is particu- 
larly important from the point of view of satisfying the 
boundary condition on the body surface, if satisfactory accu- 
racy is to be achieved on a coarse grid. 

It is possible to discretize the coupled system of 
Eqs. (3.3.4) and (3.3.5) on the physical plane, using rectan- 
gular panels with sides parallel to the co-ordinate axes. The 
drawback with this procedure is that the special care necessary 
to cluster mesh cells in the leading edge region, makes the 
grid construction airfoil dependent and complicates input data 
specification. Moreover such a mesh leads to uneconomical 
storage and Inefficient coding , especially for lifting flews. 

A curvilinear grid generated through well defined analytical 
and nimerical transformations does away with these difficulties. 
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3 . 4 . 1 . Grid Generation 

The following sequence of transformations used by 
Mercer et.al. C 102] is defined to construct the v - © plane from 
the physical x - y plane, see Figure 3.4. 


x ” r t 

y 

= y 

(3.4.1a) 

“1 

cosh ( l 

- 2 exp ( z ) ) 

(3.4.1b) 

s 

© 

= t/t 

' upper 

(3.4.1c) 

s + it 

and 

e = x + iy 



The principal transformation of course is Eq. (3.4.1b). It 
unwraps the airfoil from the trailing edge and stretches it on 
either side of the t-axis, on the s - t plane. r t is the 
leading edge radius if the airfoil is round nosed or a small 
positive value if the airfoil is sharp nosed. Eq. (3.4.1a) 
simply shifts the origin from the leading edge of the airfoil 
to the centre of curvature of the leading edge. This is 
necessary because the second transformation is singular at the 
origin which must therefore be excluded from the computational 
dcxnain. 

The last transformation (3.4.1c) yields a rectangular 
domain on which the upper and lower surfaces of the airfoil are 
stretched into 0 = 1 lines on either side of the ©-axis. The 
upstream infinity maps on to the origin of the v - © plane. The 
downstream infinities of the upper and lower surface wakes 
along the positive x-axis on the physical plane, correspond 





upper 


Fig. 3.4. Grid generation through co-ordinate transformations 














Fig. 3.5. Grid on the Physical Plane for a NACA 64A410 
Airfoil . 
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respectively to the negative and positive infinities of the 
o' -axis in Figure 3-4. 

The rectangular domain on the v - 9 plane is divided 
into smaller rectangular cells. The rectangular grid is now 
mapped back to the physical plane yielding constant v and 0 
lines/ which constitute the curvilinear grid on the latter 
plane; Figure 3.5. All derivative calculations are carried 
out on the v - 0 plane and the derivatives on the physical 
plane are obtained using the Jacobian of the transformation, 
according to the formula 

^ (3.4.2) 



where t u (s) denotes the profile surface on the s - t plane and 
prime denotes differentiation with respect to s . Since the 
transformation (3.4.1b) is analytic on the x - y plane except 
at the origin, its real and imaginary parts satisfy the Cauchy 
Riemann equations . Hence 
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at _ _ as 

3 x “ 3y 


and 


as _ at 

3x ~ 3y 


(3.4.5) 


Substituting in Eq. (3.4.4) and carrying out the matrix multi- 
plication we obtain 


a s _ 0t u ^ s ^ as i at 

ax ” t u (s) ax t u (s) ax 

3_t 6t A (s > at i a_s 

a x "t~Cs7 ax + t u (s) a x 


(3.4.6) 


The derivatives 3s/ 3x, at/ ax and t^(s) are all 
evaluated using second order accurate central difference formulae 
The Jacobian J is thus evaluated for each grid point and stored 
in memory. 

3.4.2. Discretization 

Equations (3.3.4) and (3.3.5) will now be discretized 
by replacing the integrals by sums and derivatives by differences 

Field Sources 

In each cell of Figure 3.5/ formed by the intersection 
of the constant v and 9 lines on the physical plane/ we consider 
all the flow variables to be uniform. Each cell is approximated 
by a quadrilateral and its influence co-efficients can be 
evaluated at every field point analytically, following a method 
outlined by Hess [l03] . The field velocity contributions may be 

obtained from the sum formulae 

0 

U p (x, y) = E E e kl (x/ y) 5 V1 (3.4.7) 

k 1 KX 
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kl — 

v p (x, y) = £ 2 f (x,y) g kl (3.4.8) 

where e and f denote the influence co-efficients for the field 
sources. Summation indices k and 1 respectively denote grid 
points running along v and © directions. The expressions for 
the influence co-efficients e and f are derived in Appendix E. 

Internal Sing u larities 

To carry out the discretization of the mean line/ the 
airfoil chord is first divided into three main segments as shown 
in Figure 3.6. Each segment is now further subdivided into 
smaller elements, according to the distribution given in the 
figure, by specifying the number of elements in each segment as 
input data. The elements on the mean line are then determined 
corresponding to the elements on the airfoil chord. Similarly 
the collocation points on the airfoil surface are located corr- 
esponding to the mid-points of the elements on the airfoil chord. 

On each element of the descretized mean line, the 
source strength is assumed to be uniform while the vorticity is 
assumed to have a linear variation. Following Basu [104], two 
point sources are placed on either end of the last mean line 
element nearest to the leading edge, while one point vortex is 
placed at the mid-point of the same element, see Figure 3.6. 
Denoting the internal singularity contributions to the velocity 
field by the subscript I, we have the following sum formulae 
for the velocities 

Uj-Cx, y) = S a k (x, y) 0 k + £ c^Cx, y) (3.4.9) 

k k 
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Fig. 3.6. Discretization of the Mean Line and Choice of the Coliocatio 
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v I (x / y) = s b k ( x , y) cJ k + E d k ( x , y) Y k (3.4.10) 

k k 

where the summations are carried out over all the singularity- 
elements constituting the -unknowns. The influence co-efficients 
a, b, c and d for the internal singularities are given in 
Appendix F. Note that the contributions from the point singu- 
larities have been absorbed in the vectors a and y . 

The perturbation velocities are therefore given by 

u(x, y) = u-j-Cx, y) + Up(x, y) (3.4.11) 

v(x/ y) = v I (x / y) + v p (x, y) (3.4.12) 

When the mean line is divided into n segments, we have 
(2n + 4) unknown singularity strengths to be determined. This 

'I 

implies that the tangency condition must be specified at (2n + 3) 
collocation points. Figure 3.5, yielding as many equations. The 
remaining equation necessary for a unique solution is obtained 
from the Kutta condition, which in the present method is satisfied 
by setting the vorticity strength to zero at the trailing edge. 

Inserting Eqs. (3.4.11) and (3.4.12) in Eq. (3.2.9), we 
obtain the following equations for the tangency boundary condi- 
tion. 

V I - F i,l Cx) "l = F i,l (X) COStX * SLna * F u,l Cx> “f - V F 

on Y = F . (x) (3.4.13) 

In matrix form Eq. (3.4.13) can be written as 

c ik l?o k Jfr k ,n T = tN i ] T +[h 1 ] T (3.4.14) 

= F ,\ i (x. ) cosa - sina (3.4.15) 

vJL/X X 


where 
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h i 


F u,i (x i > u f (x i / ^ - V x i' y±> 


(3.4.16) 


and (x^, y A ) denote the collocation points on the airfoil surface, 
; The vector h^ arises because of the presence of the non- 
linear term. Setting ^ s 0, gives the solution for incorrpress- 
ible flow. The matrix C remains fixed throughout the iteration 
scheme and hence may be inverted and stored at the beginning of 
the iteration. 


3.4.3. Non-linear Source Term 


From Eqs. (3.2.5), (3.2.18) and (3.2.22), the field 
source term may be expressed as 

9 = % 4a ( §” ) + T (3.4.17) 

V 

The first term in Eq. (3.4.17) represents the inviscid 
part of the field 'source distribution. The streamwise deriva- 
tives occurring in this term may be evaluated directly, but the 
following procedure for calculating them enables the governing 
equations to be satisfied explicitly. Since 


3 ,o; 
k 2 




u + v 

U ^x V Tx 


_L (£.) « u + v 

ay l 2 J u ay ay 


we obtain 


q iL (£L) 

ds { 2 } 


c r> + c r> 


(3.4.18) 


Hence if u x and u y are known, using the irrotationality condition 
and the governing equation we -Obtain 

and 


v = u 
x y 


-old 
v *= q 

y * 


u 


X 


(3.4.19) 
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The derivatives of u are computed on the v - Q plane using 
second order accurate central differences at all interior points 
and one sided differences at all boundary points. The deriva- 
tives on the physical plane are obtained using Eqs. (3.4.2) 
and (3.4.6). 

The second term in Eq. (3.4.17) represents the artifi- 
cial viscosity. Assuming that the streamlines on the physical 
plane do not deviate very much from the ^-co-ordinate lines in 
the supersonic regions for the airfoils to be considered/ we 
may use upwind differencing directly on the physical plane. 

Since the i-indexing starts from the trailing edge on the upper 
surface and increases continually through the leading edge to 
the trailing edge on the lower surface# upwind differencing 
amounts to forward differencing on the upper surface of the 
airfoil and backward differencing on the lower surface. Such a 
sinplified differencing procedure for the artificial viscosity 
term works satisfactorily for the airfoils with small camber. 

may therefore be computed 

surface 

(3.4.20) 

(3.4.21) 

is to be used, a similar 


The artificial viscosity contribution 
from 

T ij = G ij “ G i+i/j on the W er 

= G. . - G. „ . 

ij i-l ,3 

whene 

— P 3 ,c? 

G = - \i -~2 Q 

1 ? 1 Inartificial compressibility 
approximation for Eq. (3.2.24) gives 
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P — “ P i+i on the 1J PP er surface 

(3.4.22) 

= P ; ^ ,(P, - P . , . ) on the lower surface 

-i-J -1-J j-j J.— X/J 

g is then obtained from Eq. (3.2.25). 

Equations (3.4.21) and (3.4.22) often produce too little 
dissipation for strongly super- critical flows. Hence the swit- 
* ching factor M in these equations is replaced by a factor X p. 
where X is simply a numerical parameter greater than one. The 
introduction of X (to be henceforth called artificial viscosity 
parameter) provides a simple means of controlling the magnitude 
of dissipation for a given case and hence the stability of the 
iteration scheme. 

3.4.4. Iterative Procedure 

A simple direct iteration scheme is used to solve Eqs. 
(3.4.11) and (3.4.12), satisfying the boundary conditions (3.4.13) 
to obtain the unknown velocity field. For a purely subsonic case, 
the scheme represents a linear operator which, by repeatedly 
acting on a non-linear forcing term to satisfy the tangency 
boundary condition, produces quick convergence in the entire 
flow field in a few iterations. For super-critical cases, with 
the introduction of artificial viscosity, it is found that con- 
vergence can be achieved even for flow fields with fairly large 
embedded supersonic regions, without the addition of any artifi- 
cial time damping terms. These terms have however been added in 
the present scheme to the field source term g to impart stability 
and to reduce the iteration count needed for convergence of 
strongly super-critical cases. 
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/ „ \ 

I£ an estimate of g is available at the n-th cycle/ 
from Eqs. (3.4.11) and (3.4.12)/ we obtain two simple linear 
integral equations for the singularity distributions cr and y 
which are determined by satisfying the tangency condition 
(3.4.13). This involves the solution of the linear system of 
equations (3.4.14). Once a and y are determined, the field 
source term may be recalculated from the present velocity field. 
This procedure is repeated till convergence is obtained in the 
entire flow field. Since the matrix C remains fixed throughout 
the iterations, it may be inverted and stored at the beginning 
of the iterative procedure. 

The following equations describe one complete cycle of 

/ n v 

the iterative scheme. From the given g ' 


u (n) 

F. . 

= 2 

1 

ki -(n) 

i e u 

from Eq. (3.4.7) 

(3.4.23) 

(n) 

V 

F. . 

= 2 

1 

T f H -(n) 

l £ ±j g kl 

from Eq. (3.4.8) 

(3.4.24) 

■, (n.) 

h 

= F u,l 

, . (n) (n) 

(x i> v, - v, 

il il 

from Eq. (3.4.16.) 

(3.4.25) 


[£o i J <n) !y 1 } (n) ] T = [CjJ - 1 [tn k } (n) +{h )t } <n) ] T 

from Eq. (3.4.14) (3.4.26) 



(n) 

= ul 

hj 

+ “f" 1 

(3.4.27) 


= .<»> 

+ v' n > 

(3.4.28) 


I. . 

F . . 

XJ 


-(nri-1) 

g ij 

is now 

obtained from Eq. 

(3.4.17) if aritificial viscosity 


is to be used and Eq. (3.2.25) if aritificial compressibility is 
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to be used and the cycle is complete. An artificial time depen- 
dent quantity represented by the sum 

, (n) (n-1). . (n) (n-1). 

“l (u ij ' u ij " v ij > 

may now be added to as obtained above. 

It is noted that in Eqs. (3.4.25) and (3.4.26)/ the 
indices i and k run only from 1 to 2n + 4, i.e. only on the grid 
points located on the airfoil surface where the tangency condi- 
tion is satisfied. It is also emphasised that the Kutta condi- 
tion at the trailing edge is satisfied by setting y^ = 0 at that 
point for every iteration. 

The initial estimate of g for the present computations 
was to simply put g^ - 0, which in the present formulation then 
leads to the incoirpressible solution. Convergence criterion was 

f r\\ 

specified in terms of maximum change in the value of the g. . 

•**» J 

over all the grid points. In the actual calculations the iter- 


ative scheme was terminated when max I A g 1 £ 0.005. This value 

ij 

corresponds to a maximum of less than 0.5 percent change in the 
surface pressure values. 


3.5. RESULTS AND DISCUSSION 

Plane/ two dimensional/ full-potential solutions for a 
few selected airfoils, obtained from the numerical scheme 
described above will be presented here. Both sub-critical and 
super-critical results are discussed, but the emphasis will be 
on the latter. First an outline of the program code is presented, 
which is followed by the results of some numerical experiments — 
which include, a study of the effect of neglecting the non-linear 
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contribution to the lateral velocity component, the effects of 
the number and distribution of grid points, the approximation 
for the field source strength and the artificial viscosity 
parameter A . This is then followed by the presentation of 
super-critical results on symmetric and lifting airfoils. The 
discussion of these results is broadly confined to noting 
features specific to the method by comparing the present results 
with those of finite-difference or finite-element calculations 
wherever the latter are available. Lastly, some typical compa- 
risons of the FPE solutions with the TSPE solutions obtained as 
in Section 2.5.5 are shown and the chapter concludes with a 
discussion of the performance of the iterative scheme. 

3.5.1. Outline of the Program Code 

The plane, two-dimensional, full-potential algorithm 
was basically incorporated into two independent program codes 
taking into consideration the fact that the curvilinear grid 
construction is dependent only on the airfoil geometry and is 
independent of the angle of attack. The grid for a given 
airfoil is therefore generated by a program package called GRID 
which also calculated the influence co-efficients for the field 
sources and internal singularities and outputs them on the 
disk as data files. These data are then input to the main 
program FULPOT which is coded to solve the problem for arbitrary, 
subsonic free- stream Mach numbers and angles of attack. Starting 
from any given Mach number and angle of attack, a series of 
cases could be run, if desired, obtaining converged solutions 
for each case and proceeding on to the next case by changing 
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either the free-stream Mach number or the angle of attack or 
both and taking the previously converged solution as the input 
solution for the new case. 

Programs were run on the computer system DEC 1090/ 
which has a core memory of 256 K words of which only 100 K are 
available to users on the interactive time- sharing mode or on 
the batch processing mode. Programs requiring larger core may 
also be loaded, but the need to avail of virtual memory drasti- 
cally increases the central processor (CPU) time since swapping 
time is also included in the latter. Most of the present 
calculations were carried out interactively on the time sharing 
mode. 

3.5.2. Some Numerical Experiments 

(a) Effect of the Approximation Made in the Calculation of 
v-corrponent of the Perturbation Velocity 

Because of the limitation on the core memory available 
to users on our computer, it became necessary to incorporate an 
approximation in the calculation of the v-corrponent of the 
perturbation velocity. The non-linear contribution to the v- 
component of velocity represented by the second term in Eq. 
(3.4.12) was calculated only for the first level of grid points 
given by j =1. The field influence co-efficients necessary 
for the computation of this contribution was thus cut down by 
about 80 percent. As a result, for super-critical flews past 
lifting airfoils, the core requirement for a typical 47 x 5 grid 
was about 120 K words which, though still above the 100 K core 
limit set for users, was found to give reasonably speedy 
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executions considering that the use of virtual memory involved 

swapping the object program in and out of the core at the end 

of almost every iteration. The effect of this approximation on 

the surface pressure distribution was tested for a symmetric 

sub-critical case of flow past a NACA 0012 airfoil at M = 0.72/ 

00 

a = 0 and the results are shown in Figure 3.7. Surface pressure 
distributions have been plotted for three cases with the v p 

ij 

contribution to v^. in Eq. (3.4.12) being calculated for 
(a) j = 1 level only, (b) for all the j-levels, (c) for none 
of the j-levels respectively. The result from Lock [105] is 
also plotted in the same figure for comparison. It is seen that 
the approximation made above, represented by case (a) results 
in surface pressure values which practically coincide with the 
unapproximated results of case (b) and both are in close 
agreement with the results of Lock. Further it is noted that 
calculation of this contribution on the surface points is 
essential, as otherwise the error in the pressure peak is as 
much as 12 percent. 

The above experiment was repeated for lifting sub- 
critical flow past NACA 0012 airfoil at = 0.63 and a = 2° 
and the results obtained were found to be similar. Though no 
such checks were carried out for super-critical lifting cases, 
it is believed that for the relatively thin and low cambered 
airfoils for which computations were carried out, the above 
approximation does not cause serious errors in the resulting 

pressure distributions. This conclusion is also strengthened 

* 

by the good agreements achieved in some of the calculated 



Woo = 0.72 ’iOt.z 0 grid : 23x6 NACA 0012 



0.6 

0.8 


v p computed for j = 1 only 

a v F computed for all j 

v F neglected for all j 

o R.C. Lock [105] 



1.2 - 

Fig. 3.7. Effect of Neglecting Nonlinear Contribution v F to the 
v- Component of Velocity. 
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results with those obtained by finite-difference or other 
methods, as we shall see later. All the results that follow 
were obtained for the case (a) approximation as mentioned above. 

(b) Grid Effects 

Figure 3.8 shows the effect of a slight variation in 
the chord-wise distribution of the grid points on the surface 
pressure for a NLR-0.ll-0.75-0.90, quasi-elliptical airfoil 
section. It is seen that small changes in the grid point 
distribution on the chord have almost no effect on the surface 
pressure distribution. 

Figure 3.9 shows comparisons of the sub-critical 
pressure distribution on the NACA 0012 airfoil obtained by using 
(33 x 5) and (41 x 6) grids, with grid extent remaining the same. 
Lock's test data for this case are also plotted in the figure. 
Grid refinement is seen to result in higher peaks on the upper 
surface of the airfoil near the leading edge. The upper 
surface pressure distribution in the interval 0.1S xg 0.9 
remains practically unaffected by grid refinement. On the 
lower surface also, grid refinement had little effect on the 
computed pressure distribution which is seen to agree well with 
Lock's results except for a small discrepancy near the trailing 
edge. This disagreement is believed to be due to the fact that 
the results of Lock are for an airfoil which was extended 
slightly beyond the trailing edge at 1.0 so as to mate it closed, 
whereas no such extension was made at the trailing edge in the 
present calculations. 








Fig. 3.9. Comparison of Coarse and Fine Grid Results : Subcritical Surface 
Pressure Distribution on the NACA 0012 Airfoil. 
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Table 3,1 

S ur face Pressure Distribution for Sub-critical Lifting Flow 


Past NACA 0012 Airfoil 

M = 0.63 a » 2° 

OD 

X 

0 

0.0005 

0*007 

0.03 

0.075 

0.139 

0.217 

0.305 

0*40 

0.503 

0.609 

0.715 

0.821 

0.927 


Grid: 33 x 5 

C T =0.328 

Ju 

°P 

Upper 

Lower 

0.9481 

0 . 9481 

0.5417 

1.1013 

0.3578 

0.7961 

-1.0184 

0.2063 

-1 .0644 

-0.1059 

-0.9370 

-0.2373 

-0.7903 

-0.2781 

-0.6478 

-0.26.97 

-0 . 5144 

-0.2339 

-0.3921 

-0.1845 

-0.2832 

-0.1303 

—0 « 1834 

-0.0731 

-0.0773 

-0.0037 

0.0865 

0.1216 

0.3128 

0.3202 


0.99? 
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For super-critical flows, it was essential to keep 
the grid point spacing in the supersonic region approximately 
uniform in the x-direction so as to ensure that the dissipation 
term remained proportional to the artificial viscosity function 
with a fixed Ax as the constant of proportionality in this 
region, since the airfoils considered are relatively thin and 
have small camber, the artificial viscosity computed from Eq. 
(3.4.20) or (3.4.22) satisfies this requirement if the chordal 
step size Ax is uniform. The interval x^ < x < x R shown in 
Figure 3.6 denotes the anticipated extent of the supersonic 
region for a given flow. To capture reasonably sharp shocks. 

Ax was fixed at 0.05 in this interval. With this choice, a 
typical grid used for symmetric cases was 26 x 6 and for lifting 
cases, 51 x 5. On the physical plane, the latter gave a grid 
extent of -0,75 £ x £ 1.30 and -1.44 ^ y £ 1.44. For sub- 
critical cases a coarser 33 x 5 grid was generally used with the 
grid extent remaining fixed. A typical grid constructed for a 
NACA 64A410 airfoil is shown in Figure 3.5. 



In all the above calculations, the field source strength 
for all field panels, except those that are adjacent to the 
airfoil surface, i.e. j = 1 row, was computed at the centre of 
the panel. For the field panels adjacent to the airfoil surface 
two different approximations for g were tried. For scheme A, 
the conputed values of g at the collocation points were directly 
used. The results presented above are those of scheme A. For 
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scheme B, a linearly interpolated value between the values at 
j = 1 and j = 2 was used. 

Figure 3.10 shows for the quasi-elliptical airfoil 
section NLR-0.ll-0.75-0.90 a comparison of the surface pressure 
distributions obtained by thm two schemes, with Lock's sub- 
critical results. The computed pressure distributions show 
disagreement with Lock’ s results near the leading edge region. 

A similar discrepancy was noticed earlier for the NACA 0012 
airfoil. Figure 3.9. For the NLR profile. Figure 3.10 shows 
that scheme B results in better agreement in the region ahead of 
the peak although the peak velocities are considerably smaller 
than those calculated by Lock. Scheme A calculations produce 
higher peak velocities but the agreement with Lock* s results in 
the region between the leading edge and the velocity peak is 
poor. 

Different ways of evaluating the field source term 

were tried to achieve better agreement between the calculated 

values and Lock 1 s results in the region between the leading edge 

and the velocity peak, but they proved unsatisf actory . One of 

the alternatives tried was, the calculation of g in terms of 

9 ■ 

the density P and its derivatives in the form g = - Q ^ (In P ) 
but this did not qualitatively alter the results. Similarly 
different distributions of grid points near the leading edge 
region also did not affect the results. Though the reasons 
for the discrepancy could not be clearly identified, it is 
believed to be due to the following. Primarily, the constant 
source approximation made for the field panels in the leading 
edge region, where the gradients are very large, is a source 



12 ? 


- 0.8 r 



Fig. 3.10 Surface Pressure Distribution on the Quasi - Elliptical Airfoil 
NLR-0. 11-0, 75-0. 90. 
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Table 3.2 

Co-ordinates for the Quasi-Elliptical Airfoil NLR-Q.ll-0.75-0.90 


x 


y 


Curvature = d©/ds 


-1.79504 

0.00000 

49.100807 

-1.79307 

0.00918 

30.8 

-1.79046 

0.01432 

21.5 

-1.78714 

0.01920 

15.79 

-1.78303 

0.02414 

12.01 

-1.77798 

0.02928 

9.33 

-1.77184 

0.03470 

7.34 

-1.76437 

0.04048 

5.822 

-1.75528 

0.04668 

4.630 

-1.74418 

0.05340 

3.681 

-1.73052 

0.06073 

2.916 

-1.71360 

0.06873 

2.295 

-1.69241 

0.07756 

1.7810 

-1.66555 

0.08733 

1.3831 

-1.63114 

0.09815 

1.05582 

-1.58661 

0.11009 

0.79811 

-1.52848 

0.12310 

0.59896 

-1.47946 

0.13236 

0.49208 

-1.45077 

0.13718 

0.44399 

-1.39976 

0.14485 

0.37417 

-1.32987 

0.15376 

0.297747 

-1.29638 

0.15748 

0.267000 

-1.24843 

0.16228 

0.229433 

-0.94118 

0.18268 

0.120369 

-0.83289 

0.18680 

0.106563 

-0.73948 

0.18933 

0.098506 

-0.6.5319 

0.19090 

0.093190 

-0.57171 

0.19175 

0.089533 

-0.49395 

0.19200 

0.086994 

-0.41934 

0.19175 

0.085257 

-0.34754 

0.19106 

0.084117 

-0.27837 

0.18998 

0.083436 

-0.14708 

0.18684 

0.083092 

-0.05444 

0.18376 

0.083498 

0.03372 

0.18016 

0.084339 

0.09011 

0.17752 

0.085101 

0.14468 

0.17470 

0.086004 

0.19748 

0.17173 

0.087030 

0.24861 

0.16861 

0.088171 

0.29799 

0.16539 

0.089412 

0.34589 

0.16205 

0.090744 

0.39217 

0.15863 

0.092156 


Contd... 
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Table 3.2 (Continued) 


x y Curvature = d8/ds 


0.45876 

0.15335 

0.094401 

0.52204 

0.14795 

0.096768 

0.58223 

0.14244 

0.099225 

0.63943 

0.13687 

0.101737 

0.69383 

0.13126 

0.1043 

0.74555 

0.12563 

0.1068 

0.79471 

0.12002 

0.1092 

0.84152 

0.11442 

0.112 

0.90033 

0.10703 

0.116 

0.95545 

0.09974 

0.115 

1.00698 

0.09258 

0.117 

1.07833 

0.08214 

0.120 

1.14317 

0.07210 

0.118 

1.20234 

0.06249 

0.113 

1.25596. 

0.05332 

0.111 

1.35001 

0.036 56 

0.0520 

1.42858 

0.02234 

-0.070293 

1.49459 

0.01079 

-0.37512 

1.52352 

0.00617 

-0.7015 

1 . 55001 

0.00250 

-1.3626 

1.58224 

0.00000 
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Table- 3.3 


Symmetric Surface Pressure Distribution Over Quasi-Elliptical 
Airfoil NLR-G.11-0.75-0.90 


M = 0.7861 a 

OD 

X 

0 

0.0005 

0.0049 

0.0193 

0.0461 

0.0811 

0.135 

0.205 

0.275 

0.345 

0.415 

0.485 

0.555 

0.625 

0.695 

0.765 

0.830 

0.89 

0.96 

0.99 


0° Grid: 45 x 5 


1.16 

1.0516 

0.5213 

-0.0808 

-0.3685 

-0.5065 

-0.5751 

-0.5891 

-0.5734 

-0*5494 

-0.5182 

-0.4797 

-0.4339 

-0.3789 

-0.3105 

-0.2233 

-0.1140 

0.0350 

0.2480 

0.4752 
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of inaccuracies. Secondly,, the internal singularity model 
adopted in the present method appears to be increasingly sensi- 
tive to the choice of 6, the leading edge gap, as the airfoil 
becomes relatively more blunt. This is illustrated in Figure 

3.11 which displays the computed results for the quasi-elliptical 

airfoil section NLR-O. 1025-0. 675-1. 375 at M = 0.7557/ a = 0° 

ao 

for two different values of 5 = 0.1 r t and 0.2 r t * This is a 
shock free super-critical case. The airfoil has a very blunt 
leading edge with r fc = 0.166 5. The figure shows that the 
surface pressure distribution is erratic near the leading edge, 
changing rather drastically in this region as 6 is increased 
from 0.1 r fc to 0.2 r^. The pressure curve in the interval 
0.2 S x S 1.0 is quite smooth though the computed values 
remain sensitive to 6 till about 0.4 of the airfoil chord. 

Thus the model gives inaccurate solutions for airfoils with 

relatively large nose radii, and actually failswhen the nose 

* 

radius is of 0(0.1) or more. Thus the disagreement between 

the computed results and Lock* s results is somewhat greater for 

the NLR-0.11-0.75-Q.90 airfoil (r, = 0.0204) because it is 

t: 

relatively more blunt than the NACA 0012 airfoil (r fc = 0.0156). 

Results from schemes A and B were also compared with 
the results of finite-difference calculations for a super- 
critical lifting case. These comparisons are shown in Figures 

3.12 and 3.22. It is clear that scheme A results show much 
better over all agreement with the conservative finite-difference 
results than scheme B. For this reason, the results to be 
presented were all calculated using only scheme A. 



! 

i 




Table 3.4 


Co-ordinates for the Quasi-Elliptical Airfoil 
NLR-0 .1025-0.67 50-1 . 37 50 


x 


y 


Curvature = d©/ds 


- 1.69058 

0.00000 

5.962550 

- 1.68599 

0.03897 

5.95 

- 1.68141 

0.05473 

5.94 

- 1.67683 

0.06656 

5.94 

- 1.67225 

0 .07630 

5.94 

- 1 . 6676)7 

0.08469 

5.941 

- 1.66810 

0.09208 

5.943 

- 1.65851 

0.09873 

5.944 

- 1.65391 

0.10477 

5.940 

- 1.64928 

0.11033 

5.929 

- 1 . 6:4460 

0.11550 

5. 904 

- 1.63984 

0.12034 

5.858 

- 1.634 97 

0.12492 

5.787 

- 1.62993 

0.12929 

5.671 

- 1.62464 

0.13355 

5.509 

- 1.61905 

0.13771 

5.282 

- 1.61286 

0.14195 

4.991 

- 1.60840 

0.14481 

4.760 

- 1.60605 

0.14625 

4.633 

- 1.59830 

0.15072 

4.225 

- 1.58935 

0.15542 

3.797 

- 1.57912 

0.16027 

3.402 

- 1.56790 

0.16505 

3.1064 

- 1.55664 

0.16935 

2.992 

- 1.54660 

0.17279 

3.111 

- 1 . 5379 a 

0.17546 

3.130 

- 1.52414 

0.17926 

1.6930 

- 1.48295 

0.18874 

0.75233 

- 1.40363 

0.20334 

0.44926 

- 1.31607 

0.21597 

0.32544 

- 1.26589 

0.22200 

0.284110 

- 1.21107 

0.22773 

0.251461 

- 1.15159 

0.23310 

0.225478 

- 1.08781 

0.23797 

0.204587 

- 1.02005 

0.24223 

0.187383 

- 0.94843 

0.24579 

0.172493 

- 0.89843 

0.24774 

0.163222 

- 0.84606 

0.24935 

0.154169 

- 0.79100 

0.25059 

0.145239 

- 0.73255 

0.25141 

0.136536 

- 0.66998 

0.25178 

0.128305 

- 0.60268 

0.25162 

0.120864 

- 0 . 54889 ) 

0.25110 

0.115972 


Contd 
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Table 3.4 (Continued) 

x y Curvature = d©/ds 


-0.49222 

0.25018 

0.111760 

-0.43282 

0.24884 

0.108267 

-0.37098 

0.24703 

0.105494 

-0.28542 

0.24387 

0.102864 

-0.19766. 

0.23984 

0.101343 

-0.10880 

0.23496 

0.100786 

-0.04222 

0.23077 

0.100911 

0.02383 

0.22618 

0.101426 

0.08893 

0.22122 

0.102308 

0.15259 

0.21594 

0.103478 

0.21513 

0.21034 

0.104903 

0.27584 

0.20451 

0.106546 

0.33476 

0.19847 

0.108370 

0.39183 

0.19225 

0.110343 

0.44699 

0.18590 

0.112434 

0.50021 

0.17944 

0.114613 

0.55152 

0.17289 

0.116848 

0.60092 

0.16629 

0.1192 

0.66387 

0.15744 

0.1226 

0.72358 

0.14859 

0.1245 

0.76037 

0.13974 

0.130 

0.83405 

0.13097 

0.129 

0.88490 

0.12230 

0.132 

0.93307 

0.11377 

0.133 

1.00051 

0.10125 

0.134 

1.06262 

0.08914 

0.131 

1.11930 

0.07762 

0.125 

1.17250 

0.06644 

0.1121 

1.22095 

0.05598 

0.0910 

1.26560 

0.04614 

0.057322 

1.34448 

0.02855 

- 0.07781 

1.41119 

0.01412 

-0.42141 

1.44054 

0.00830 

-0.7959 

1.46746 

0.00362 

-1.5593 

1.50134 

0.00000 







Fig. 3.12 Surface Pressure Distribution on a NACA 0012 Airfoil , 
Scheme B . 
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(d) Artificial Viscosity Parameter 

The artificial viscosity parameter X is significant 
from the point of view of convergence of the iteration scheme 
because, once the artificial viscosity form is chosen, the 
parameter X provides a simple means of controlling the magni- 
tude of the dissipation term in the supersonic region at super- 
critical Mach numbers. 

Figure 3.13 shows the effect of the viscosity parameter 
X on the surface pressure distribution calculated for symmetric 
flow past a NACA 0012 airfoil at = 0.82. For this case, X 
was initially set at 3.0 and convergence was achieved in 50 
iterations. The calculations were repeated with X = 2.6 and 
2.4. For X = 2.4, the iterative scheme showed signs of going 
into an oscillation and was terminated after 100 iterations; 
the surface pressure distribution at this stage is plotted in 
Figure 3.13. 

Figure 3.13 shows that the peak velocities attained just 
ahead of the shock increase as the value of X decreases. The 
shock location remains more or less fixed at about 0.575 chord. 
Converged solutions could still be obtained with X =2.5, but 
with a further decrease in X to 2.4, instabilities occurred. 

The surface pressure plot shows the incipient formation of a 
peaky region ahead of the shock. The formation of such a peaky 
region however does not necessarily lead to divergence . But it 
is a sure sign of the occurrence of instabilities which may or 
may not arrplify. Further decrease in the value of X however 
leads to quick divergence with the appearance of a prominent and 
growing over expansion region just ahead of the shock. 



a = 0° Mocr 0.82 grid: 26x6 


X 

— 2.6 

2.4 (Unconverged) 


Effect of X on the Surface Pressure Distribution for a 
NACA 0012 Airfoil . 
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Fig. 3.14. Effect of A on the Surface Mach Number Distribution lor a 
NACA 0012 Airfoil. 
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This is clearly illustrated in Figure 3.14 in which 
the surface Mach number distributions are plotted for a slightly 
higher free- stream Mach number 0.83. The converged solutions 
were obtained for X = 4.0 and 3.0. The surface Mach number 
distribution obtained at the end of 20 iterations for A =2.8 
shows a sharp and prominent peak in the shock region which 
finally led to divergence of the scheme in about 35 iterations. 

Hence it appears that solutions differing significantly 
in the shock region, but agreeing broadly elsewhere, are 
obtainable depending upon the value chosen for the artificial 
viscosity parameter A . Such arbitrariness appears to be 
associated with other in viscid, potential flow formulations 
too and is the result of the need to introduce artificial 
viscosity in the numerical schemes, see for example [36]. Some 
circumspection is therefore necessary in choosing the appro- 
priate solutions. In the present thesis the following criterion 
has been adopted: From among a sequence of solutions for 
decreasing values of X , that solution is accepted which does 
not show a peaky velocity profile ahead of the shock and which 
corresponds to the least value of X * Therefore, it becomes 
necessary to generate solutions for atleast two or three values 
of X starting from some appropriate high value (depending upon 
the airfoil, free- stream Mach number and angle of attack) and 
decreasing in steps varying from 0.2 to 0.5* In the present 
work this has been done for almost every case and surface pres- 
sure distributions were obtained for atleast two values of X • 
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3.5.3. Symmetric Super-Critical Results 

Symmetric super-critical solutions with shocks were 
obtained for flow past a sharp nosed parabolic arc airfoil and 
the NACA 0012 airfoil at different free stream Mach numbers. 

Figure 3.15 shows the surface pressure distribution 
obtained on a 10 percent thick parabolic arc airfoil at = 

0.84. Solutions obtained with X = 4.0, 3.0 and 2.5 are 
compared with the results of conservative finite-difference 
calculations by Holst [l3]. The present method is seen to 
predict higher velocity peaks as compared to Holst* f results/ 
but the shock locations predicted by both methods agree well. 

The finite-difference shock is seen to be weaker than that 
obtained from the IEM. 

Figures 3.16 and 3.17 show super-critical symmetric 

solutions for flow over a NACA 0012 airfoil at free-stream 

Mach numbers ranging from 0.78 to 0.84. Figure 3.16 shows the 

surf ace pressure plots while Figure 3.17 shows the corresponding 

Mach number plots. The X value quoted for each free-stream 

Mach number is approximately the value at which converged 

solutions, satisfying the criterion cited in Section 3. 5. 2d, 

could be obtained. The Mach number plot shows that shocks are 

captured better for higher free-stream Mach numbers, these being 

sharper and better centred, as judged from the pressure plots 

of Figure 3.16, than the shocks at lower free-stream Mach 

numbers. Thus for instance, the shock at M = 0.78 is seen to 

oo 

be practically smeared out. 

A comparison of the present results with those of 
conservative, finite-element calculations for two cases. 
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Fig. 3.15. Surface Pressure Distribution for a Parabolic 
Arc Airfoil. 





Fig. 3.17 
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Fig. 3.16. Comparison of Surfact Prossuro Distribution *ith Fiilftp- 
Eiomont Conssrvativs Rt suits for a NACA 0012 Airfoil. 
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Fig. 3.19. Comparison of the Surface Pressure distribution with Finite? 
Element Conservative results for 0 NACA 0013 Airfoil. 
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Table 3.5 

Surface Pressure Distribution for Super-critical Symmetric 
Plow Past NACA 0012 Airfoil 


M = 0.85 a = 0° X = 4,0 a- = 2.0 Grid: 26 x 6 
ao 1 


x 



0 

0.0005 

0.0041 

0.0193 

0.0461 

0.0811 

0.125 

0.175 

0.225 

0.275 

0.325 

0.375 

0.425 

0.475 

0.525 

0.575 

0.6.25 

0.675 

0.725 

0.775 

0.83 

0.89 

0.95 

0.91 


1.1939 

1.1078 

0.6352 

0.0173 

-0.3458 

-0.4997 

-0.5981 

-0.6751 

-0.7264 

-0.76)53 

-0.7950 

-0.8185 

-0.8356 

-0.8463 

-0.8512 

-0.8486 

-0.8367 

-0.8086 

-0.7529 

-0.5027 

0.0113 

0.3281 

0.3994 

0.5227 
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illustrates the shock capturing characteristics of the present 
method. Figures 3.18 and 3.19 show comparisons of the surface 
pressure distributions obtained by the present method with the 
results of Bristeau [37] for = 0.80 and 0.85 respectively. 
Figure 3.18 shows that the captured shock for the present method 
is rather smeared as compared to the shock obtained from the 
conservative, finite-element calculations. Further the shock 
location obtained by the present method is about two mesh widths 
upstream of the shock of [ 37 ]. Figure 3.19 on the other hand 
shows that the present shock location agrees fairly well with 
the one calculated by FEM although it is stronger than the 
latter. Similar trends were observed for the lifting cases and 
these results will be discussed in the next section. 

Lastly it is remarked that in Figures 3.16 and 3.17, 
the converged results for M ^ = 0.84 could be obtained only by 
including the artificial time-dependent term a-jU^ to the field 
source term g with = 0.5. The effect of a 1 on convergence 
history is discussed in a later section. 

3.5.4. Super-Critical Results for Lifting Flows 

Super-critical solutions for lifting flows were obtained 
for the symmetric NACA 0012 airfoil and the cambered NACA 64A410 
airfoil whose profile is presented in terms of the upper and 
lower surface co-ordinates [109]. 

Figures 3.20 to 3.22 show the surface pressure distri- 
bution obtained for the NACA 0012 airfoil at angles of attack 

= 1.0, 1.5 and 2.0 degrees. Values of A for which the surface 
pressure distributions were obtained together with the calculated 
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Fig. 3.22. Surface Pressure Distribution on a NACA 0012 Airfoil 
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Table 3.6, 



Surface Pressure Distribution for Super- critical Lifting 
Flow Past NACA 0012 Airfoil 


.75 cc = 2° X “ 

3,5 cc^ = 0.5 

Grid s 51 x 

554 C M = -0,1515 

Opj — —0.0118 



c 


X 

p 



Upper 

Lower 

0.0 

0,9319 

0.9319 

0,0005 

0.9276. 

1.1486 

0.0049 

-0.1361 

0.9146 

0.0193 

-0,7713 

0.4092 

0.0461 

-0.9633 

0.0579 

0.0811 

-1,1299. 

-0.1216 

0.125 

-1.2140 

-0.2297 

0.175 

-1.2708 

-0.2886 

0.225 

-1.2978 

—0.3104 

0.275 

-1.3109 

-0.3111 

0.325 

-1.3156 

-0.2964 

0.375 

-1.3128 

-0.2728 

0.425 

-1,3032 

-0.2433 

0.475 

-1.26-54 

-0.2129 

0.525 

-1.1320 

-0.1783 

0.575 

-0.4677 

-0.1499 

0.625 

0.0686 

-0.1193 

0.&7.5 

0.0117 

-0.0808 

0.725 

0.0100 

—0.0648 

0.775 

0.0044 

-0.0176) 

0.845 

0.0458 

0.0414 

0.935 

0.1772 

0.1794 

0.99 

0.3845 

0.3843 
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C L and C M values are also Indicated in the figures . The moment 
co-efficients were calculated about the origin. The + sign in 
front of a number under the column ‘iters' / indicates that the 
iterations for the corresponding value of X was continued 
from the converged solutions for the previous value of X . 

Figure 3.22 shows that the computed surface pressure distribution 
for M = 0.75, a = 2 ° compares rather well with the conservative 
results obtained by finite-difference calculations [99] parti- 
cularly as regards shock strength and location. However for 
M = 0.75 at a lower angle of attack, a = 1.5°, Figure 3.21 
shows that the computed shock is considerably weaker and is 
situated about one mesh width upstream of the shock calculated 
by another conservative finite-difference scheme [ 106] . Thus 
here again, as for the symmetric cases considered earlier, the 
present scheme seems to give results which show better agreement 
with finite-difference conservative results as the flow becomes 
more super-critical. 

Figures 3.23 and 3.24 show the results obtained for a 
NACA 64A410 airfoil. Figure 3.23 shows the surface pressure 
distributions at three different free- stream Mach numbers, 

0.735, 0.74 and 0.75, while Figure 3.24 shows the plots of the 
corresponding surface Mach number distributions. The rapid 
increase in the extent of the supersonic region and the 
strength of the shock as the free- stream Mach number is incre- 
ased from 0.735 to 0.75 is clearly evident. The Mach number 
plot for M^ = 0.72 is also included in Figure 3.24. This is a 
super-critical case and Jameson’s finite-difference, quasi- 
linear calculations [25] show the presence of a sharpshock on 




Moo 

X 

C L 

Cm 

0.735 

2.0 

0.617 

-0.29 

0.74 

2.5 

0.643 

-0.30 

0.75 

3.0 

0.865 

-0.424 


Surface Pressure Distribution over a NACA 64A410 Airfoi 
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Table 3 0 7 

Co-ordinates for the NACA 64A410 Airfoil 


Upper surface 


x 


y 


0. ‘JO 90 2 


0 

0.0035 
0.00582 
0.01059 
0.02276 
0.04749 
0.07230 
0.09737 
0.1474a 
0.1977 
0.24a 
0.29834 
0.34871 
0.39910 
0.44950 
0.49989 
0.5502 5 
0 .00057 
0.65085 
0.70108. 
0.75126 
0.80151 
0.85148 
0.90104 
0*95053 
1.0 


0,01112 
0.01451 
0.02095 
0*03034 
0.37&5 
0,04380 
0,05366 
0.06126 
0.06705 
0.07131 
0.07414 
0.07552 
0.07552 
0.07344 
0.07040 
0.06624 
0. 06106> 
0.05490 
0.04780 
0.03957 
0.03018 
0.02038 
0.01028 
0,00021 


Lower 

x 

0 

0.00650 

0.00918 

0,01441 

0.02724 

0.05251 

0.07770 

0.10263 

0.15252 

0.20230 

0.252 

0.30166 

0.35129 

0.40090 

0.45050 

0.50011 

0.54975 

0.59943 

0.64915 

0.69892 

0.74874 

0.79849 

0.84852 

0.89896. 

0.94947 

1.0 


y 

o 

-0.00678 
-0.00796 
-0.00969 
-0.01251 
-0.01592 
-0.01819 
-0.01996 
-0.02244 
-0 .02406 
-0.02499 
-0.02537 
-0.-2518 
-0.02436 
-0.02266 
-0.02024 
-0.01736 
-0.01418 
-0.01086 
-0.0076 
-0.0046 
-0.00229 
-0.00132 
-0.00076 
-0.00048 
- 0.00021 


surface 
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Table 3,8 

Surf ace Pressure Distribution for Super-critical Flow Past 
NACA 64A410 Airfoil 


M = 0.75 

CD 

& 

If 

O 

o 

X * 3 

. 5 Grid : 

47 x 5 

P 

11 

o 

* 

C L = 0.845 

C M " 

-0.413 

C D 

-0.016 






C 






P 


X 



Upper 


Lower 

0 



1.1206, 


1.1206. 

0.0005 



1.0512 


0.9088 

0.0049 



0.2581 


0.2153 

0.0193 



-0,1467 


0.0817 

0.0461 



-0.3947 


0.0006 

0.0811 



-0.5712 


0.005 

0.125 



-0.6791 


-0.0114 

0.175 



-0,7712 


-0.0266 

0.225 



-0.8489. 


-0.0386 

0.275 



-0.9185 


-0.0477 

0.325 



-0.9910 


-0.0534 

0.375 



-1.0738 


-0.0588 

0.425 



—1.1486> 


-0.0511 

0.475 



-1.2032 


-0.0159 

0.525 



-1.2443 


0.0255 

0.575 



-1.2720 


0.0675 

0.625 



-1.2880 


0.1114 

0.675 



-1.2863 


0.1553 

0.725 



-1.2352 


0.1997 

0.775 



-0.9688 


0.2419 

0.83 



-0.2799 


0,2407 

o.sa. 



0.1887 


0.2221 

0.95 



0.2275 


0.2608 

0.99 



0.3384 


0.3456 
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shock on a 256 x 64 grid and a smeared shock on a 64 x 16 grid 
at 0.5 chord. The present calculation on a very coarse mesh 
of 47 x 5 has practically smeared out the shock. 

The last example of a super-critical calculation. 

Figure 3.25, shows a comparison of the results obtained for a 
NACA 0012 airfoil, M = 0.75, a = 2 6 using artificial compre- 
ssibility and artificial viscosity. The use of artificial 
compressibility as given in Section 3.2.3 is seen to produce a 
surface pressure distribution (C L = 0.45) which is closer to 
the non-conservative results (C L = 0.444) than conservative 
results (C, =0.58). 

Li 

The CPU time needed for the NACA 64A410 calculations 
was about 3 seconds per iteration. For the NACA 0012 profile, 
it was about 1.5 seconds per iteration. The computational 
effort however was the same for both cases. The difference was 
found to arise because of the use of virtual core for the NACA 
64A410 calculations which caused swapping time to be added to 
the execution time, resulting in larger CPU time per iteration. 

The maximum local Mach number reached in these examples 
is slightly more than 1.3 which is the upper limit of validity 
of the weak shock assumption. Getting the iteration scheme to 
converge for still stronger super-critical flows was found to 
be difficult. The iteration count increases rapidly as the 
flow becomes more and more super-critical. There are however 
reasons to believe that with the systematic use of artificial 
time-dependent damping terms, it may be possible to attain 
converged solutions for super-critical cases stronger than 
those presented here.- . 
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3.5.5* Iterative Development of the Solutions and Convergence 

Performance of the Scheme 

Figures 3.26 and 3.27 illustrate the convergence perfor- 
mance of the iterative scheme for symmetric flow past a NACA 0012 
airfoil . 

Figure 3.26 shows the iterative development of the 

surface Mach number distribution at M =0.84. The value of ^ 

ao 

was set at 4.0 for this case. Artificial time-dependent term 
a-j_ u t had to be used to achieve convergence in this case, with 
cc ^ = 0.5. For smaller values of o^, the iteration scheme tended 
to diverge. For = 0.5, convergence occurred in about 110 
iterations. The figure shows the surface Mach number distribu- 
tion plotted for every 20 iterations. It is seen that the shock 
location has stabilised itself around 0.70 chord in about 20 
iterations and thereafter changes very little. At the end of 
about 40 iterations, an over expansion tendency is noticeable 
just ahead of the shock. This however does no damage to the 
further development of the solution which, to plotting accuracy, 
has reached its converged value in about 60 iterations. The 
rest of the iterations are carried out to achieve finer conver- 
gence based on the criterion on A 5 * see Section 3.4.4, 

The plot of A g for this case is also shown in Figure 
3.26. It is seen that the convergence behaviour for this case 
is not uniform and that there are large fluctuations in the 
initial stages. These fluctuations die out in about 50 itera- 
tions after which, convergence though slow is fairly smooth, 
j; 

Figure 3.27 shows the effect of the addition of an 

u 

artificial time-dependent term ttjU^. to the field source term g 
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Fig. 326(conw). Iterative Development : Mach Number Distribution on the 
Surface of a NACA 0012 Airfoil. 
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Fig. 3.27. Variation of Number of Iterations for Convergence with 
(a 2 =0) for a NACA 0012 Airfoil. 
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Fig. 3.28. Iterative Development : Surface Pressure Distribution on o 
NACA 0012 Airfoil. 
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during the iterative calculations, with a 1 = 0, the iter- 
ative scheme diverges for this Mach number. Convergence occurs 
in about 110 iterations when =0.5 and as is increased 
to 1.25/ convergence is achieved in just about 50 iterations. 
With further increase in ct^ however, the iteration count begins 
to increase. 

Figure 3.28 shows the iterative development of the 
surface pressure distribution for a lifting flow, the airfoil 
again being NACA 0012. For this case, M = 0.75, a = 1.5° and 
= 0.5. Here again we notice a tendency for peaky velocities 
to develop ahead of the shock at the end of 20 iterations. The 
shock position has not stabilised yet. The solution at the end 
of 40 iterations when Ag £ 0.025, is however indistinguish- 

able from the final converged solution after 49 iterations. The 
lower surface pressure distribution has expectedly converged in 
less than 20 iterations. 

The plot of Ag shows that the convergence rate has 
declined as the solution is approached for A g < 0.5. The 
convergence on the other hand for this case is more or less 
smooth. This is due to the moderate size of the supersonic 
region and the relatively weak shock. 

3.5.6. Comparison of FP and TSP Solutions 

A few comparisons of the FP and TSP solutions obtained 
for a 10 percent thick parabolic arc airfoil will be presented 
in this section. Both FP and TSP calculations were carried out 
without the use of artificial time-dependent terms. 





120 iterations 
9? CPU (sees.) 


Surface Pressure Distribution on a Parabolic 
Arc Airfoil. 
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Figures 3.29a and b shows the computed FP and TSP 
pressure distributions at two different free- stream Mach numbers, 
0.825 and 0.85. The agreement between the two results is seen 
to be rather good. The FP shock is located slightly downstream 
of the TSP shock and is also somewhat stronger. Figure 3.30 
shows the converged FP solution obtained for M = 0.875 whereas 

QD 

the TSP code failed to converge for this case, even with the 
inclusion of artificial time-dependent terms. 

The computing times for the FP and TSP solutions for 
M = 0,825 appear to be comparable . However it must be kept 
in mind that the FP code, for symmetric cases, uses about 4 
times the core required by a TSP code and hence is four times 
as costly. This advantage of the TSP code seems to be off set 
for M ^ = 0.85, where it takes nearly 4 times as much CPU time 
as the FP code which converges rather quickly. This, as well 
as the fact that converged solution could be obtained for 

= 0.875, suggest that the FP formulation, for the IEM, is 
numerically more stable than the TSP formulation. 



CHAPTER IV 


CONCLUSIONS 


From the results and discussion presented in Chapters 
II and III, it is possible to draw the following conclusions as 
regards the use of the Integral Equation Method for the compu- 
tation of axisymmetric and plane, two-dimensional super-critical 
flows with weak shock discontinuities. 

In general, for both TSPE and FPE formulations the 
shock-capturing field panel technique of IEM seems to work 
satisfactorily for flows in which the embedded supersonic regions 
are large and terminated by relatively strong shocks. The 
captured shocks in such cases are sharp and agree fairly well 
with comparable results of FDM or FEM, in terms of both shock 
strength and location. For flows which are mildly super-critical 
the field panel technique tends to smear the shock, and both the 
strength and location of the calculated shocks differ consider- 
ably from the FDM or FEM results. The situation here can 
perhaps be inproved by resorting to a finer mesh but this has 
the undesirable consequence of increasing the core requirement . 

The wall interference calculations for the axi symmetric 
(TSPE) case show that the IEM is also well suited for dealing 
with different boundary conditions, e.g. perforated wall, slotted 
wall, solid wall, free jet or specification of a known pressure 
or velocity distribution at the tunnel wall in the problem. 

These boundary conditions then appear as additional integrals 
in the integral equation for the axial velocity component. 
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Modifications to the program code may therefore be easily 
incorporated in both TSP and FP numerical algorithms, to take 
into account different wall conditions. 


Convergence rates achieved for the present straight- 
forward algorithms may be considered satisfactory. Even for the 
strongest super-critical cases for which computations were 
carried out, convergence was achieved in about 120 iterations. 
Considering that a detailed numerical study of the convergence 
characteristics of the present schema has not been carried out, 
it is reasonable to believe that improvements in convergence 
rates are possible If more attention is paid to this aspect. 

For the TSP version of the axisymrreatric problem, two 
iteration schemes the quasi-Newton Schema and the Direct 
iteration Scheme - were used. Results of the present investi- 
gation suggest that for free air calculations, the QNS is 

faster and more stable. For the. T.,-,n ^ 

* me wall interference case, the 

Dis appears to be more stable than the QNS. 

Unlike the semi-direct methods which use fast Poisson 
solvers and generally diverge for even slightly super _. critical 

flows without the use of artificial time-dependent danping 
terms, the present integral equation approach, especially the 
FP version for plane flow, has shown that the inclusion of 
artificial viscosity is sufficient to achieve convergence for 
most cases arising in inviscid, isotropic flow. However, the 
inclusion of artificial time-dependent terms has been found to 
be beneficial for achieving stability and convergence for 
strong super-critical cases. This is true for both the 
and FP formulations presented in the thesis. 


TSP 
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A theoretical analysis of the effect of the presence 
of the time-dependent terms on iterative convergence has not 
been attempted in the present work. Therefore the range of 
values of the parameters co , and <v within which the itera- 
tion scheme converges is not known in advance. This restricts 
the usefulness of these parameters, at present, to only obta- 
ining convergence for strong super-critical cases which other- 
wise might diverge . 

The viscosity parameter X in FP formulation is found 
to be significant for its effect on the strength of the captured 
shock. The criterion proposed in the present work, by which a 
solution is chosen from a sequence of solutions obtained for 
decreasing values of X , is to be considered tentative , since 
it has not been tested for a wide range of airfoils. The problem 
here, as already noted, is that the artificial viscosity term, 
being a completely arbitrary numerical artifice, there is no 
satisfactory consideration by which the amount of dissipation 
can be prescribed in advance. In the absence of suitable 
physical criteria, this term is manipulated only by numerical 
consideration such as stability or captured shock characteris- 
tics. Hence, small and at times significant, uncertainties 
are present on the computed surface pressure distributions. 
However it must be noted that this uncertainty is not a pecul- 
iarity of the IEM alone. It is common to other potential flow 
calculations also, sea for example [36]. 

In the present approach to the FP formulation, the 
internal singularity model of Basu [l04] for incompressible flow 
has been successfully adopted for the calculation of compressible 
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flows with or without shocks, by considering the non-linear 
terms as an additional field source distribution. The advantage 
of the internal singularity model as compared to the surface 
singularity models, particularly for purposes of compressible 
flow calculations with shocks, lies in the fact that the former 
yields satisfactory 1 ncor.pressible solutions with fewer number 
of unknowns. Clear!/ tl:rr, the singularity model used to 
generate the incompressible solution to the Laplace equation 
must be satisfactory if the final non-linear solution is to be 
satisfactory. As mentioned in Chapter III, the present internal 
singularity model fails to yield satisfactory results for 
airfoils with nose radius of the order of 0.1 or greater. 

Further the calculated solution in such cases are sensitive to 
the choice of the leading edge gap 6. For these reasons it 
is clear that there is a need to either modify the present 
model or search for alternative surface singularity models 
which retain the advantage of the internal singularity model of 
requiring fewer unknowns, see for example [110] . 

.Lastly, it must be recognized that, IEM algorithms 
generally need, for a given problem and mesh size, several times 
more core than a corresponding finite-difference algorithm. The 
difference in the core requirements arises largely due to the 
need to store the influence co-efficients. The TSP and the FP 
versions of the IEM algorithms, respectively need for a given 
grid, 2 to 3 times and 4 to 6 times the core required by the 
corresponding finite-difference algorithms,, Some core economy 
can be achieved for the case of TSP algorithms as conpared to 
11? algorithms because the former calculations are carried out 
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on rectangular grids which give the influence co-efficients 
certain symmetry properties. Even so, for a given problem and 
mesh size, the EDM may seem attractive purely from the point of 
view of core economy. However some of the features of the IEM 
enable an overall reduction in the total number of grid points 
to be achieved, which reduces the core requirement considerably. 
Firstly, as was evident from the discussion of both TSP and FP 
results, the IEM gives satisfactory results with a much smaller 
grid extent for the infinite domain problem than the finite- 
difference and finite- volume methods. Secondly, to realise the 
same accuracy, the IEM computations can be carried out on a 
rather coarse mesh and therefore need fewer grid points than 
the FEMs. For these reasons, the number of grid points used in 
IEM are generally 1/2 or 1/4 the number used in equivalent 
finite-difference or finite- volume methods. 
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APPENDIX A 


EVALUATION OF EXPRESSIONS OCCURRING IN SECTION 2.3 


(i) Sinplif ication of tha Expression for I 


sw 


From Section 2.3, subsection (iiib) we have 


sw 


- r w n cjzr v. +1 * d ©‘ d * 


p A *x ! r 1 =r 


(A-l) 




w 


00 2 % 


w 


P -CD 0 


/ ^ L?| d®‘ dx* - 


r* =r 


w 


CD 2% 

r w ' ' s' Sr- 1 

- CD 0 


d©» dx* 


r 1 =r 


w 


Now consider each integral separately over 0* with ¥ ** 


and R = [ (x - x* ) 2 + r 2 + r 2 - 2rr w cos6‘] . 

1 F(ft/2, k) 


4n R 


Thus we obtain. 


2% l % dfi 

' * de ' - h ' * 


It r / t \ 2 " , I > 2-] 1/2 

[ (x - x' ) + (r + r w ) ] 


(A-2) 


4rr, 


where k » 


w 


(x - x ' ) 2 + (r + r ^) 2 


(see 291.00 of ref. 107 ) 
Again 
2Xi 


/ s , de* 

v r* 


, 2% r - r cos©' 

•±= / -A « d0* 

4,1 0 R 3 

% w 

■*- [ (r - r) / + r / = d©'] 

w 0 R 0 R J 



183 


1 r (r~ _ 2E (71/2/ k) 

271 L w 77 .'. 2 , ' 2-, l/2 r " . .2 , “27 

L (x-x 1 ) +(r+r w ) j[ (x-x‘ ) +(r-r w ) ] 


w 


, 1 ,P(V2, k) - E(n/2, M -. 

r w [(x-x' r + (r+r w ) j 1//2 

(see 291.01 and 291.06 of ref. [107] ) . 


(A— 3 ) 


Hence the expression (A-l) leads to 

r - r 


QD 


/ E(r. W 


E(7t/2, k) 


w 


- QD 


W % 


[(x-x«) 2 + (r -r) j [ (x-x 1 ) 2 + (r+r 

W W 


2X Ux-x') 2 + Cr+rJ 2 ] 172 


■5T 7 *T ,/2< i, T71 <*x< (x ' ' r w )] dx ' <^ 4 > 

” [(x-x') 2 + (r+r..) 2 ! 1/2 x 


w 


(ii) Differentiation of the Wall Term in Eq. (2.3.10) 2ith 
Respect to x 
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APPENDIX 3 


INFLUENCE CO-EFFICIENTS FOR AXISYMMETRIC FLOW 
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(i) Derivation of the Field Influence Coefficient a_ ( x,r ) 
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where F is the complete elliptic integral of the first kind and 
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The elliptic integral is now treated as a constant 
and taken out of the integral sign with its value at the centre 
of the element. Therefore,, 
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(ii) Derivation of the Wall Influence Co-efficient a (x # r) 


From Eq. (2.4.4) 
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The wall co-efficients consist of a symmetric point 
arising from the terms A and B and an anti-symmetric part 
arising from D„ Treating the elliptic integrals as constant 
over each of the small elementary intervals of integration, the 
resulting simple integrals may be easily evaluated. Denoting 
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the symmetric part by a g and anti- symmetric part by a A the 
integrations are carried out to obtain 
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expression for a g becomes singular. It is therefore separately 
evaluated as follows. 
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~~ % f>o 6 + ~ 3 ~~ ^ + ^ + 2(ln 6* - 1) S b Q 6 + 6 


+ — 6 * 5 } ] 
D 


where 6 = 6/V 


Hence for x = x -, , and r = r w we may write 


a k (x,r) 


B 2n 


! ^ + [X 2 + (r + r w ) 2 ] 1/2 

E(k) i n - ^ * 2 t l/2 


Xi + [ X x + (r + r w ) ] 


The co-efficients a Q , etc. and b Q , etc. are [108] 
1.3862944 a, = 0.1119723 a« = 0.0725296 


0.5 


a x = 0.1119723 
b x « 0.1213478 


b 2 = 0,028879, 


The wall influence co-efficient is obtained from 


a w (x,r) « a s (x,r) + a A (x,r) 


(B-5) 


(B— 6} 


(B-7 ) 


Par a cotrpu'- ; v ; „na! grid with uniform mesh width in 
the axial direction, tJv* following relations of symmetry hold. 


kl 

F U 


il 


F 


kj 


XI k-1,1 

a. 23 ■ 

ju' ' 

*44 * 4^*1 4 
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Similarly 

k i , 

a„ = a_ ana 

S . . Si . 



a*' 1 and 




a*' 1 
A. . . 

3 .- 1 / J 


Hence the storage ar-“ .-required for a p is (k x 1 x j) location 
and for it is (2 x k x j) locations. 



APPENDIX C 


EQUATIONS FOR SLENDER BODY GEOMETRY 


The equation for the profile of a general parabolic 
body of revolution is given by the following two equations, 
ref . [ 96] . 

For body with maximum thickness forward of its mid-point 
R(x) = C R mrx [1 - x - (1 - x) n ] (C-l) 

For body with maximum thickness aft of its mid-point 

R(x) = C R max (x - x n ) (02) 

The following table gives the data for various C and 
n combinations with locations of maximum thickness. 


x at R 

max 

C 

n 

Formula used 

Sting 

0.3 

1.71 

6.03 

C-l 

0.71 

0.4 

2.36 

3.39 

C-l 

0.78 

0.5 

4.00 

2.00 

C-l or C-2 

0.85 

0.6 

2.36 

3.39 

C-2 

0.89 

0.7 

1.71 

6.03 

C-2 

0.93 


The cross-sectional area function s(x) and its derivatives are 
s(x) = it R 2 (x) 
s* (x) *= 2 it RR* (x) 

s" <x) * 2 % (RR u (x) + R ,2 (x>) 






APPENDIX D 


CALCULATION OF ARTIFICIAL VISCOSITY FOR PLANE FLOW 


Simplif ication of artificial viscosity term T, Eq. 


(3.2.18) 


T = - Jg [ As v Q (f" ) ] 

c 


'Q 


(3 .2.18) 


o r . <- . a . 

•g-g LAS Ji — g t U (' 


.3 +V*> +v 4 (£+£),] 


3y 


3 L A SU™{(U 2 U +uvv)+ (v 2 v + uv u )}] 


■Si 


1 a , „ 3 v r Q p #.,2. 


13 (u + v a? 5 [ TJ Ax n “7 < u u v + uv 


“2 

C 


X 


+ ~ Ay p. (V 2 v + uv u i )] 


(D-l) 


Inspection of Eq. (D-l) shows that its highest deri- 
vative terms are the same as those in Eq. (3.2.15). Neglecting 
the cross derivative terms in (D-l) we obtain 


3G , 3h 
3x + 3y 


(D-2) 


■where G 


A X U -4 K + UV V J 

c 


and H » - Ay P< (V 2 V + VU U ) 

c y y 



APPENDIX E 


CALCULATION OF FIELD INFLUENCE CO-EFFICIENTS FOR PLANE FLOW 

Influence Co-efficients for a Plane Field Source Element 
of Unit Strength 

Figure (Evl) shows an arbitrary, plane, quadrilateral 
source element S of unit strength. The co-ordinates are 
referred to axes fixed at the grid point of the element whose 
corners are numbered 1, 2, 3 and 4, going clockwise. The 
co-ordinates of this grid point in the global co-ordinate 
system fixed at the leading edge of the airfoil are ( rj-^) 

The u and v components of velocity, in the x and y direction 
respectively, induced by S at an arbitrary field point (x, y) 
are given by 

ff * dx‘ dy* 

S (x-x* ) Z + (y-y‘) Z 

(E— 1) 

// Y . ,Z Y . ! dx* dy* 

s Cx-x«r + (y-y*) 

(E— 2) 

The integrals may be conveniently evaluated by writing 
each as a sum of integrals over 4 infinite strips. First we 
consider the evaluation of the integral in Eq. (E-l) . 

Consider the side 1-2 of the quadrilateral element 
1-2- 3-4-1 shown in Figure (E-l). This side is now assumed to 
be replaced by a source sheet of infinite extent with sides 
parallel to the x-axis as shown in Figure (E-2) . Assuming 


u g (x,y) = e^ (x,y) = 

v s (x,y) = f^x^y) = 
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that the quadrilateral is traversed in the clockwise direction 
as indicated in Figure (2-1), we assign a uniform source 
strength -1/2 for the semi- inf inite sheet lying to the left of 
side 1-2 and +1/2 for the semi-infinite sheet lying to its 
right . 


Traversing the quadrilateral clockwise/ the other 
sides are also replaced by similar infinite source strips with 
source strengths +1/2 or -1/2 according as the semi-infinite 
part of the strip lies to the right or left of the side being 
traversed. Thus one infinite strip is associated with each 
side of the quadrilateral and the source strength undergoes a 
jump of unit strength across the side. These source strips 
may be denoted by where 

v — p, + 1 for (i = 1/2/3 

= 1 for p, = 4 

If denotes the source strength function having 
values +1/2 and -1/2 on either side of the line p - v, we may 
write the integral in Eq. (E-l) as 


kl , 

e (x,y) 


4- £ V 

2n 


{1=1 (x- 


X - X* 

:-x’)^ + (y-y* ) 


* — 5 a \xV ^x' dy« 


; kl, . 
2 QnV (X/Y> 
R=1 


(E-3) 


Each integral in Eq. (F/-3) may be explicitly evalu- 
ated. For example/ consider the integral over the strip S^y. 
The equation for the line joining |i- v is given by 



{ fr-4 ) 
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so that 


x* = p + q , y 
F )XV \iv J 


where 




\iv 


i 

V u 

jj. 


- <v 4 + 4 

kl 


The integral for e'~ can be written as 


'pv 

y‘ 


3tL , s 

V (x ' Y) 


v 


4t 


/ dy‘ J - -i / 


i *• <y > i 


CD 


(2-5) 


( 2 - 6 ) 

(3-7) 


y 1 

P 


- CD 


+ "9 / 

2 x‘ (y* ) 


} 


IF 


In 


[ (x— x l ) 2 + (y-y‘) 2 l dx‘ dy* (s-8) 

where x’ (y‘ ) is given by Eq. (E-5). Eq. (3-8) holds for both 
Yjl > Yy and < Yy • 

The integral in Eq. (3-8) can be easily evaluated in 
terms of elementary functions and the following expression may 
be obtained 

ejj£(x,y> = Xk [(y i " y i Un h 3 + F(y *' Y li' h 2' ^ lE ' 9> 

where 

F(y, 6 , a, P) = (y + a) m [(y +ot ) 2 + p*i - 

(6 + a) In [( 6 + a ) 2 + p 2 ] - 2(y - $). + 

-1 y +a -l 6 + a % ir 

2P (tan - tan — 

Dropping the subscripts pv t the quantities h^ h 2 ' h 4 


axe obtained from 
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= C (x - p ) 2 + y 2 ] /hi, 
h 2 = (x - p) q + y A 13 

o 2 

hf = 1 + q 


^2 *2 
h l h 2 


■E 11) 
(H-12) 
(£-13) 
(E-14) 


The calculation of the influence co-efficients for 
the v-component proceeds in the same fashion except that in 
this case, the strips have sides parallel to the y-axis. 
Replacing y‘ by x 1 , x by y and vice versa we get. 


-kl , v 

V (X/y) 


[(4 - x‘)ln g 2 + F(^, x', g 2 , g 4 >] (2-15) 


f^"(x,y) = 2 fj^(x,y) 

p=l ^ 


(B-16) 


The corresponding expressions for g^j g 2 etc. are as follows. 


y. - 



X* - 

X* 

(E-17) 

y ^ — 

J v 



'K - 

X 1 

U 


Y l 

=5 

2T 

\xv 

+ 

X 1 

(l>18) 

"C2 

S3 

y 1 

— y 1 


(E-19) 

S \lV 



' x n 




S3 


+ yi 

E-20) 

2 

g l 

S3 

r 2 
[x 

+ <y - zrt/vl 

(.0-21) 

g 2 

«3 ' 

[x ■ 

4- (y - 

t)s] /g^ 

(3-22) 

2 

g 3 


1 4- 

2 

s 


(3-23) 
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2 2 2 
g 4 “ g l “ g 2 


(E— 24 ) 


The function F Is the same as defined by Eq. (E-10). 

Co-efficients calculated using Eqs. (E-l) to (E-24), 
though exact/ involve large number of function evaluations. 
When the plane source element is far from the field point it 
can be approximated by a point source and the influence co- 
efficients can be evaluated from the formulae 


e^Cx/y) 


1 _ 

2 % 


x - x* 


(x - x 1 ) + (y - y* ) 




(E-25) 


f kl U / y) 


2 % 


Y ~ V 1 


(x - x' ) + (y - y‘ ) 


1 A 3cl 


(E— 26) 


where represents the area of the quadrilateral element 

k,l. Formulae (E-25) and (E-26) were used when the condition 


[(x - x* ) 2 + (y - y') 2 ] 1/2 > 4 

was satisfied, with denoting the length of the longer of 
the two diagonals of the source element. In actual calcula- 
tions it was found that about 40 percent of the total number 
of influence co-efficients could be calculated from Eqs. (E-25) 
and (E-26). 

The co-efficients e and f calculated as above remain 
the sane for the airfoil fixed co-ordinate system since the 
induced velocities u and v are invariant with respect to siiqple 

translation of the origin. 

To economise on storage for special cases, the foll- 
owing symmetries in the influence co-efficients may be noted. 
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(a) Flow Past a Symmetric Airfoil at Zero Angle of Attack 

In this case for every plane source element S on the 

★ 

upper plane/ there corresponds a conjugate element S on the 

★ 

lower plane. The geometry of S is the mirror reflection of 
S. Because of symmetry, the calculation for such a flow needs 
to be carried out only on one half of the plane. Hence the 

-k 

velocity components induced by S at a field point on the 
upper half plane must also be taken into account. Since the 
induced velocities for two source elements conjugate to each 
other are symmetric in u and anti- symmetric in v, we have, 

u *<x, y) <= u s (x, -y) (E-27) 

s 

v *(x, y) = - v (x, -y) (E-28) 

S 

so that the total induced velocities at (x, y) are 

u(x, y) = Uj, (x, y) + u *(x, -y) (E-29) 

S 

v(x, y) = v s (x, y) - v *(x, -y) (E-30) 

s 

(la) Lifting Flow Past Symmetric Airfoils 

When the airfoil is symmetric, the computational 

grid on the physical plane is symmetric about the x-axis. It 

is then obvious that given a source element S on the upper 

★ ■ 

half plane, there corresponds to it a conjugate element S on 
the lower half plane. The velocity conponents induced at a 
point (x, y) by S* of unit strength will therefore satisfy 
Eqs. (E-27) and (E-28). This implies then that w© need to 
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store only the influence co-efficients of the source elements 
on the upper half plane. These must however be evaluated for 
all the grid points both on the upper and lower half planes. 

For general flow about an arbitrary cambered airfoil, 
neither of these symmetries hold. Hence for this case, the 
entire array of influence co-efficients have to be computed 
and stored. 



APPENDIX F 


INFLUENCE CO-EFFICIENTS FOR THE INTERNAL SOURCE AND 
VORTICITY DISTRIBUTIONS 


(i) Uniform Source Distribution 

Referring to Figure (F-l ) , we consider a source 
element k of width 6 and choose co-ordinate axes with origin 
at the mid-point of the element and x-axis coinciding with the 
element. Bars below refer to the element fixed co-ordinate 
system. 


The velocity components induced at a field point 
(x, y) by a source element of unit strength are 


u(x, y) 

. +6 T 

1 r x - 5 

d I 

(F-l) 

~p *. 2 -2 

-6 (x - § ) + y 

_ 1_ + ; Y 

v(x v y) 

dF 

(P- 2 ) 

2% -l liDfTTF 

Carrying out the integrations we 

obtain. 


u(x, y) 

, ,- ..2 -2 

a 1 in S 2 L±j£+JL- 8 
4 lT r c v 2 -2 

(x - o) + y 

= F(x, y # 6) 

(F— 3 ) 

v(x# y) 

1 1 x + 6 -1 x - 6 \ _ 

« ( tan — Z — “ tan — 3 -— } *=* 

■O' (x, y, 6) 


If the element inclination to the x-axis of the air- 
foil fixed co-ordinate system is 8 and the co-ordinates of 
its mid-point are (x Q , y Q ) , we obtain the influence co-effi- 
cients in the global co-ordinate system as follows. 

(x - x^) cgs 8 + (y - y Q ) sin8 




(P- 5 ) 
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Fig. F-1. Influence co-efficient calculation for the line elements 
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y = - (x - x ) sine + (y - Y 0 ) COSQ 


a^(x^ y) = F (x/ y, 6)cos0 - G(x, y, 6)sin© 
b k (x/ y) = F(x/ y, 6)sin0 + G(x/ y, 5)cos0 


(F-6) 

(F-7) 

(F-8) 


where x and y airs co-ordinates of the field point in the air- 
foil fixed (global) co-ordinate system. 

(ii) Linear Vorticity Distribution 

A linear variation of the vorticity strength is 
assumed over the element k in terms of the values at its end 
points . 

Again assuming an element fixed co-ordinate system/ 


Y k+1 + Y k Y k+1 " Y k „ 

Y (x) = 2 23 


(F-9) 


The velocity components induced 

given by 



u(x/ y) ** 

1 *r 

Ip* 

y - 

251 - 6 

(x - 1 ) 2 + : 

v(x/ y) =* 

l 

X- F 

' ^ .6 

"(x - !)" + 


Y (T) df 


(F) dl 


(F-10) 


(F— 11 ) 


Substituting Eq. (F-9) in (F-10) and (F-ll) carrying out the 
integrations and transforming to the global co-ordinate system/ 
the influence co-efficients c and d may be obtained as (n 
below refers to the number of mean line elements) 


c k (x/ y) - A k-1 (x, y) + B k (x, y) 


(F-12) 



203 


d k (x/ y) 

= C k ' 1 (x^ y) + D k (x/ y) 

(F-13) 

c^’Cx, y) 

= B 1 (x, y) 

(F-14) 

n+1 , . 

c (x, y) 

= A n (x, y) 

(F-15) 

d 1 (x # y) 

= D 1 (x / y) 

(F-16) 

d (x, y) 

= C n (x, y) 

(F-17) 


A (x, y) = -g[(G + - •^)cos0 + (P — -g— - ^e— + -^j.) sin©] 

(P-18) 

B k (x/ y) *= -^(G ” + ^)cos9 + (F + ^ - ~) sin©] 

(P-19) 

C k (x, y) = 1 [(g + ^ - ^)sln0 - (f - |£ - |2 + i-)cos©] 

(F— 20) 

D k (x / y) • |[(G - ^ +^)sin6 - (F + ^ ^ - ~)cosQ] 

tF-21) 

where (x/ y) refer to the co-ordinates of the field point 
while F and G are functions defined by Eqs. (F-3) and (F-4) 
with arguments given by Eqs. (F-5) and (F-6). 


(iii ) Point Sources 

The velocities induced by the point sources of unit 

v It 

strength may be included in the co-efficients a and b . 
we have 
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b^(x/ 

(iv) 

at (x, 

j 

c k (x, 

dNx, 


y) = 


y - Yi 


2 % , ,2 , , ,2 

(x - x k ) + (y - y k ) 


k = nrf-1. n+2 


(F-23) 


Point Vortex 


The induced velocities due to a point vortex located 

k k 

,* y^) may be included in the co-efficients c and d . 


1 Y " *k. 

y) = Tte 7 — 2 

(x - x k r + (y - y k ; 


3 k.. = n+2 


x - 


y) = 


*k 


_ y " """ T2 

(x - x k ) + (y - y^.) 


(P-24) 


(F— 25) 



i {.sms 


■SF fUg PR 35 k* 

■ T ' ' . 


vn*0- rw 

m 




r 

C 

c 

r 

C 


r 

r 

r 


r 

r 

r 

c 

r 

r 

ft 

c 

c 

c 

c 


V. 

r 

C 

c 

c 

r 

r 


c 


Mi) i s t r vcasirn n r r.^r ^ns^A 1 ** 

l-r ]f>v lA^til r \r urmS TSV' T'J JMJ i)Vt£ SAMPLE SXECU F T 3« .* 

” -i ? *, o r o o r an s to" 1 '1 he executed t o o ? f. o o r with 
ago.- d.p'nc, u GdD.'d”, AHP^'rp.FiR an j A^r/s 1 .FJR 
31.1 al si i Lied I Li hag UHnrv, 3efnr a executin'! 
r i« s ^.nivi tie 1 1 f 1 1 a w r od vnn *ii c coefficients “’SH and r S u 
Tint "M.wtM 0” “K-ecitl hi 1RID.F1R with *-«F;SS.KhR,A B :R e ’3Ii. r 3R 
•’i»5f> caef.fi dents are outbjt in t.hg Us*? 4 n the files 
" D d n „ rv, ” if) i r :j R 2 1 . n A T , fie na in nrn ?ra 1 ' is now r »ady 
f a r » x n c 1 1 1 1 a n , 

ills )T]ir j n c n m n j t e s llftlm compress* lhie flaw aver airfoils. 

~ )ivir"»eic n if the !*■ erat inn scheme mav he expected for maxlmj«r 
i ara 1 dd numbers un to 1.4. tie Imir data format as well as 
rid 1 * n e a n i m 1 s exalalneJ in tie lata flip R 1 « n a r . 'is er 
iVeri"t.hn incomes necessarv wien the it erati nn llilt set ov 
T r •■ h v is exceed?*. Anar poor late messaies are gl vei at such 
1 1 t e s ail ns°r na v r A snpni according to nis choice. 

* r cnn^ei-i a nce stiller ncssaies ir* disol ave 1 nn the terminal 
an ' us »r r»smn is wi r i anor par t ate Hants. 

tVs in ■* \ T u pfprtran th 3 t Includes subroutl nes UVFSIG, HVFSF2 

rs'-CAL, hFRraP.hyPlR and «* IRC*.- fhe rest of tne subroutl nes 
r ° c .i 4 re ) ci in f r n n the routines Tentlonea aop ve , The output Is 
written ok me piste In the fli» c '3R2i .DAT, The subroutine 
s' j i A a F Is a 'J Ah rantin'* use! for matrix Inversion. 

dp subroutine Add?,. P 3R has Its own lata file for inout.Tni s 
lira fid np-f,.nAT contains lata pertaining to the airfoil. 

o A "> S 0 A '* nevEf.oPCh M n V4 MOATS’! d : W.S.RAVirHANDRAN 

to tie oest of the authors cnodeloe the program Is free of 
ojos tor tie airfoils for wnfch calculations were carried oat 
r amort r 1 in mis thesis. However me program has not been 
r nour on ihi v chpcve^ for all tne opt ions .in particular the 
s jhro jtioe AFRS*pr ( ias not been t*st run for all. possible ways 
1 1 ohin the aerial i coordinates may oe prescribed, Therefore 
f p r s pme of the cases not tried so far ougs will very 1 lltel y 
b2 tojni in this subroutine. 


hATEjJ'HF 10 , 1 HR 1 


h m EH 51 1H VS!?(55) .TS'K'j?! 

D I U E *1 5 T 3 U WKAWFAf IhO) 
h l HE v ST 3 f 3 VCH («551 .YCHC^SlfCPCHCSI) 
hiHEvSTIU ZETAfS?) 

?3vpr,gx T, GRID, TUFT 

DIHSNST3H GP'ICSS) .^PM *551, CPl» I iy(5§) 
DI*»EHST3H CPMM(*>51 ,CP(55) 

DIWEMST3H rcf5D;sa>;CCIMtf{*53,59) 


I m •. m mw •* m — m l)f *»«*«* « ww *» <* <* m n* *• # #•* 


;PLtNL(55) 






c 

r 

c 

fm 


-;}vir i / \> / j? { ;<i 1 , ,')*X (51 *3) , ^¥(51 ,5) 

* vj 1 1 >i /'*?'/ ~ ‘ i V i 1 { 5 t , 5 ) 

ijif'm'. ■* <* j, 

5‘U U i J v< 


•> ?* i 


v? ;f i'l [ Pal , nt-;'';t^*'^K%FTL? = r B1 .OATM 

9 « ' If i UV" 

>> K i > 1 4 U, JC,SFJ,X IbAST^irTMTH.DITLV, 

1 < Srjr,, vsn^ , , SP5TE 

•? •-■ > r i , MS«S, I m, ALP1 ,»LP2 


■J U r P ( "> ; , 7 1 1 1 
? )}*.*, r *■/ * . y ■» it i'hi/i 

p n nrflM ' Vi-*,n (MhJMST DV*T. T*AMS 3*ir PAST t» L R N 61 URFnir.s'/ 
i J- r, r r!i,r..p-m:'irnr 1 amd f. * a r t tmgemcy nonvnwv 

■>2 1 j) T P T f) '! S )V TRP TJTrTRH 5TUATT j'l ' E T HDD , ' / // / 

r 0 i 1'idh'i *s the input 3ata so*»ci f IcstionM 

nf ■? T r P C 0 J , 7 9 1 1 

If < T fF(“> ),f>n2)A7HC, 

T )? *A ff I S< ' AtfPJa' f t.D , „ 

4 < r r F c ^ 3 . n ; ' n fsm , alf a , cav , lamax 

«■?* V&rr/f'»V'rsHa > ,PJi,S, 3K '' AbF Aa% F6. 3, 3X, 'GAMS' r F6.1,1X, '1,4*1 AX* 

4 # |J* ^ ,*j| yf ^ 

)3S5Ii??tf 5^^h^;'5^ 5 In,S X ,'S F1 =',r5.2,5*,' N ^ST.-,T3/! 

s*3P 4AT( >X , 'hSTUT tts^.FS. 3,*U>h£LX*'»F7.4 ,5X » 'nELTAa # ,FlO„7/) 


P3P 4 A ft I'iX. 'Fr’imSs ,Fn, 1,S1 ,'jELXs ,f ’.‘IpHf I't'il** ,1* iw»* 

1 >6.3/) 

»RrrF(?)»701) 

calculate the basic constants »n^ parameters of the flow. 


tmmmm-m m *» » ■»' » » p **> ; 


n tP*l!t 41 1 p2?51s 5^7P12? EXPSA*«1 • * 2l , -?* P ? 9 t f/PIR 

SSSS^SKSiS-’l.-Ejift L, * i ^ ;|f J j” 

75bMTC»S / )R T C t l r )/A'>jf A fj FAD* ALFA) J ALFA* ALJAf GAMI *l. + 5*0 
"PSrARsPRF*!! 1 .-F3AM*tOSOMTC**2«t *))*»( JEXP.*AJ)-1 .0) 

4L^aALFA*t»i»/!90.» £SALFA*C3 S( ALFA) ) SnALFA*STN( ALFA) 


C set up the nr Is* 



?JSSt*.8?iSS?8!:5ftSE«*attleXa:SiAi:«aK:Sa;»i:KE?KdlllR B 

compute the coefficients for the derivative celculetion 

?|i^ 

compote tne f nf iutnce coefficient* for the internel •in^uietlti* 

- -i ■:, h 1 ? a ■: - 1 .* . r ••'•: . n.,'i J , « n* & , *« \ < » 7 Jlas r,l^) 

tit j s ^ ^ n / } 




,< x r p ( ? ) ,'i'W '* i ’ai.'J, fJ « U" , ’b, v a* 

.» j T ,' 5 -f •» ) siro , t\y -< iw,^, mx ,x‘’r!,v^K 

..• ,r ?)/,o ^ >5 mr,r«ri'R 

,i ; f 1 1 'i 

^ m ^ m ' mmm mm mm 

n \ f , aTl Mp’ri 

t :>" , f ri|T=^, ^ur: I '=':> r >s*) 

^ ? • , r i • » 1 1 * s ? i, \ y ? v r : r = ' •) «; k ' ) 

Vi i:> ?=i ,*mx 

S?r. >} ijo, iois f (*S'l(T f 

irvlf ?rx m '\lH f ( 2 sV(T,*,f,) »K=l ,<'*R v ),li=l,i’ , A») 

*■*■) v | ,| i ? r 

j i r r k ( r > ,' * i *; ° ) 1 

- j ' 1 rr rrt; 

s- •] a a v ( \ r > r i i , » ) 

r ■)’> m t f ? y , " i ntmt ’-ni ’* e a :l f'r 1 - . l“) 

r ■) i> i o 4 n h 

i ■>=• j r j’t i Ts-> \ , i^’tCPs'D"'^ ' ) 

?SA ^ (ASfK,!, I*] | ' « j 1 5 2 

>f 2 J*> )'(*$£<, I. Jj.jsj » PAll’ s 'a ,W 1 

5 SMf ) f ("v < , I , J f S , M J< , S , JX 

cj £ft Df 2 ? , «v p ) r { r* V f < * I . J) , Ts 1 r * » Tst * M 4R X 1 

nv'‘ 

F'JSMA r fim 3. 7 1 
"D'UTVli „ ^ 

~Mu RTl'fSf 1*^3) _ _ _ 


'alruiat* <■">'* "a! locfltbn matrix 


ID 1 
D J 1 
TFf l , 

". 2 ( 1 , 

22 ( 1 , 

2 d n 
2 2 ( 1 , 

22 ( 1 , 
mi. 

T ? f l , 

civr! 

r 3 *j r t 


3 T*1 f }IC 

??i«gfK^}*^In- 5 fPCCn<'A 5 fKrl + 'JWr 1 ) 
ID? 

n*AS(K. I + !«W,1 ) 

X f •iljS + ?)xC'l(K, I + MW. 1 ) 

£?.' M CD C?C T *1 . Kj=0,0 
£?>C)C?CT'M f K*ML^?)a«0.0 


xns 

M'JS 


?(Tft f MLS + 1 )*l, 


t avert the natrlx CC(I ,tO 

2\hl ' raiAAFCCC*! A#Nm!cCI«V, IA, ^«CARCA,inST) 

." ^ ■ ■ m ^ m ^ jig. i,.' m m m m mm# >ng up -in ''fill m wd -Mtt •* "• •* in ii *ji <■* • W w 1 

rbA 3 «,»At*SR. 

f >3 IS Tsi i sklX .. ; ? . ! ; * 

d s «irv.i , , 




'4 V c T , << ) * /,r [ h £ < ) 
t ?) 'MVMftf l-^fjS-2 


jf« r i« ► a ► a l velocity f ^ 2 ? i 


MC'M®y= J” J”?* 

'•^r.ruT )-\ .- : !C“! i t> 2-vrjMp? 

*’ ) i r t <!i>: 

f T P I a 1 * f T K f » s M '* I j(MsO g O 

- vil, fTR'I,TrP( 1 ,X'',v5,C°liTN f \T ) FIV,YH,Cr lf rD,C^) 

mm mmmmmmmmmmmmmmmm m m m mmm m m m mmmmmmmmmmmrnmmmm 

r !•! ( ^ s , i n 

•j u r k ( ? :> , 5 u •» c r - , r [> , x x , r * 

<^TrE(^a t 31 3) 

VJ 1 1 Ts’ , C 

•nr IT (2 J, sr,M <"( T ) , C n LiT t ' 

" 3’H’T NM£ 

-R* rr(?>,7<',i> 


rrai. cu late tr>° M° l* sources 
"■} *£ r t M 1 T £ 

%r,u 3F! S T V C U p , ?FX , UFY ) 

Vj ?t 

') 3 72 T*t , MM S x 

’3 ** X ( T » 1 ) a U C i f l . J ) 

mUsrs*S' f, U6' p CT, l 1) ? '(/CEJMpaSSftLr^iVrU* J) 
T(T. J)sjrcWP**? + v Cn f «P**2 

^H;H(T, 1)3S'JRT(fKT,.J)M2(T, 1)) 

ilnf f 5i2ri?o2£*6rxcx, jnvco*p*mct, J) 

Jj)-uft)U*EXPcl« STER ” 

iffiiiisifeW 1 *"--'-- 




T , 1 J s;i«* If r , I 1 + 0„05*(.GF'IC T, 2)-G p, if l, H ) 

j 11 i' T 4 ’ 

3 " f T 2 
'JV1*T V'f;: 

'i o a *•’ s n . >: 

1 i = i X 

IT, T)sr,nsrf;p-if ! , 

=• r r i , n r , jn 

i , , ’) 

*i’? rr v?*; 

1 j ’! FT '4’?2 

k i e *i t — 1 *>* =/ ft v • 

’ ? f n =P RF * r ( 1 . ~fn A“ * r 3 n , i *» * * 2 - 1 . y ) ) * * f 3 » **EXPG a r ) - 1 . 0 l 

■ r '( 4 o '3 r n , 1 i).' j. o iRrr&r 2 ^ , R !j M T r . C r H , T s 1 ^MftX) 

nr TfivsviMT.BM^M. K it T r F C s f *1 h ) T P , R N 5 R 

? ip* A If] *X, ‘TTs*, n,rX,'R FIS. 7) 

rrr r*i i/r cp<;}n j rn 1 v n 1 J 1 F ( R 'JOR v .HT. ^0 . 0 ) ,J m 10f>2 

risiTfl; Tf rir.CT.r’MftXloin n l^bJ 3"> T 0 107 


!IJ5fI7*m$u» trttrH^.n.oun «jlom L< 

lonpute tne induced velocities due to $ne fiell so-irces and 
i> *if>n calculate t n e non* 1 1 neat f otci fi9 function vector -F. 

:nLb JWFSF7 
J ft L L F FGCAT, 

B ' ? '*cnb calculates trie f orclno function at toe collocation points 
TFf IT.fV tlCftLt, RlTMF£tSF34) 

C 3 rn u>3 

oasic Iterative sequent ends 


;JmFc5375?2nf7ft^RM, Sn rn ia53 

^TrF( 29 , 57 mt f RNnRM 

"pmpj 

"?( n «**{♦} ? ? ?-FSA#* £3(1# D **2-1 ■,- 0 >)-**CCM*eKP 6 * , «) - 1,01 

*3*ir intis 

03 03 TatjMMAX ,»% 

xsiiCT>«*FAr.t?G5l2(T £ i}> 


FS'JCDs^Fftr.lFG^iniT.U J 

V 5 UT)sftTnftGfZGRT 3 f 1 , 1 )) 

nnurrMUs 





a n a a 


•ui>:>»rs'Vi'U» of attflcc/viseositv oaraneter . 

fioru’nn ^ara^tsrs nr.pl , « 1^2 -nay also be varied lenendina 
:ij«^ r -> a jiprf of tairp or toe i 1 ft r a 1 1 v * scheme so far. 


c 

r 


JRTJ FIS 5SR) 

3 '• a f 5 , * > I " r. r. 

T «’ f I " i\ r 1 , ^ J « u 1 S T 3 3 

sriur s.-nmy.Af f r< 
^n^S^USbPi.M.p^r.ft^x 

PiVjr s. o nr.pi ,M> 2 ,LM*.nx 


" nr- L, =?T 

oiosl.1 
PX&'iMs 
Nfs?,/ 
■) S *1 't T " ~ 
" o «; r A 3 s 

•jutfc? 

^ t r s - 1 0 
,nT r so ( ”> 
*’nx r c? 
v-»rr r (.? 


^1 502^515097912; EVPj^ = 1 • / £2?!!* , A , 2i r-ft 0 1 *PC»0 
. / r - R *! - 1 . ) ; F 5 M 2 sFS m ** 2 » FCU«sn . 5 * f SAM -1 ,0 *Fs«? 

S M / p R m ? • A 3 *D„ 5 *(i . ♦ 3 *M 1 ; n 3 * t , / p 5 w 2 + 0 .5 * C .. A M -1 .1 

tny r ( r\Q / & 0 ) * M*F ft D»M*F ^ # T IU» F & = ALFA • GAM 1 st* + Gft®l 

?r, + " 1 «»/? 0 n. j CRftf'FRsC35f ftLFftl ; SNM,rA»STNf a£f*) 

0 , 5 f 0 1 

*) r S031FSh f T»U, M-F»D 
L»,bf'9)050»'ir,rpSTM? 

Dl552-ur.P1 , R(,P2,liAM&K 

’•> , 5 f n _ _ __ - 


I r* AV=5, 



fSBiM 


s-DOaR.rf 2 V , 1 ? 0 c 1 H-) ) 

*39M»rnrx;'*cm'.8X,'CPLTNCI)'l 

s , 3P5»jf i«:x,F10,4 r 5X,Fl0,4) , „ c i/itv '.wenli 

if /I 5^, 'free s f r e « t. t*acb nu^orr FSM- »f5»^J5X#alrroil 
rb-i»cF 0 PSS T A” = ',r5,3/i5X f '’anpie of at tac£ ALFAfi* ,^5.2/ 

USX *arMf icm viscosity factor f _ 'I 6 - 21, „ 

rr-DRMa H / 1 5v , 'number of line elements NL.S- , T3'l ox , 
l 'na-noor or collocation joints wr* ,131 

U55i5?; J Sh5r° f i"ti?5ai‘5sML< rS'SKS 

PDPlJn/15X f 'son1c ve|ocltv 3s’)MTC*',F7.3/15X # 'son!c pressure 
s-DR^Aff 27 iTMSx^iimi t^exceedfi , tyoe new ITMAX or tyoe sero 
P 3 RiAr f I5x! 'uncooveraetl results J r *! » n»HAl5S2# M r« 

;353:ffli5^*SSSH*2«5SaiU5 r {l5.*?|6l5fe;Jf97 s -tSlfsSi.-* 

:3!«f ;?k :iih;s; 1 ; ? ; lW¥*i'wnu» 

s? ssh sr.«iawm:*l?5.TK; / « - 

/i5X?'tvp» new nu»g«r»n'i angle ot r atta^ U 

HS5s;?5|?;?icii^TitiJS sssssfitu**^? • »•«>«/ 

h«s dll 



******* * ******* *M***** ******** ************* 

CAL 

“ ■ . A- ~ .tijiirtijiiii’ .! . - , * 1 1 •, 


***************** * 
t; J n ^ «VJ T I V |- ij vr p *? f f? 


********************** 


Tnl s s n or on tine calculates tee induce. 1 velocities due to toe 
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S j'« v r BiJ »; </ + ns f K , t , J 1 *.*5 J 0 ( K } + o V ( K , T , J ) * 3 A M A ( K ) 

C 3 u J* T u u r 

HFPCT, 1) =IJF( T,.1) J VCPCI, JlsVpa, Jl 

'.if a , j ) ssu»»u+uF«r> c t , j j ? f ( t , j) =sumv ♦ vfnl n , Jl 

C3«ri¥UE 
C3¥ rTNMS. 

R £ T (J R M 
FMO 

*************************************************************** 
S J n r n j t i M K UVFCF? 

Tnls snorontlne calculates tbe non-linear contribution to 
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tapoo 1 of the t ransformat Ion. 

o I v £ k? s t o v o(«5 1 ,51 r ox( 5 1 1 5 ) a 0 T C 5 1 

C j 'J v; n-ij /vi / K S.NMAX, M A * # * If J » 
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SUBROUTINE-:, DERCOF 

This suorout tne computes the coefficients for derivative 
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SUPR0UT1ME FHRCFCH , T2 , X , Y , CP, ALPHA , Xd , CD , CD , CM ) 


compotes lift- drft»7 and moment coefficients by trapezoidal, rule 

CM=/>.0? Naf2-l 
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FUNCTION FPUX) 
FPLaO.O? RETURN 
END 

***************** 
FUNCTION FPtit(X) 
FPLl* 0 , 0 ? RETURN 
END • 

****************** 
FUNCTION FPCCX) 
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RETURN 
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<un<n.n , rie s » l* r f r ^ r n „ s, f,pp, * p? f pppp^^vm, k»j, vn, midp;, t no ) 

tragic s^Une 

sol ine is fitted to data array F at. nodes s fro® Index * to 

fad 3 vr 0 

<Hsl , 2, i-indf cates that f irst second or third derivative is 
o i v e n v i i a e v 3 a t o o i n t m 

= l , 2 , 3 — 1 n.-i •* cates that first second or third derivative is 
oiven v a 1 ue v.'j at ooirit <1 , . 

M '■moF=0 nodal values of first second and third derivatives 
of; soilne are stored in F-P r FPP»FPP» arrays so that fitted 
value at a distance H beyond a node is 

?' t .df'cpfd**?*r , pD/?. + f^*i*Fpop/d,. 

ff nin > o , fppp is given the nodal values of the integral of 
^'instead of its third derivative, starting with the value FSd 
tnen the third derivative can ne recovered as 
f f’9?f I+n-FPf I))/(SCT + 1 3 "SCO) 

r d n is set enuai to zero if s is not a monotone array 


hlMS'iSIdd sc 1.50) ,FCi5i)),FPCl50) .FPPC150) ,FPPP(150) 
r'n = h; k*i n«iS(M-*ij ,• ifck-i 5*1, 31 i 

K = f l - d) /K; T = M; .T = M + K 

hS=3f J)-S(T)? h=DS 

r Ff dS) 1 1 , 8 1 ,11 

h F 3 C p C 1 ) " F f O ) /OS 

'1 = 0, 5; v = 3 . £ ( 0 P - V tf ) / b S ’ 33 ff) 25 

'1 = 0,0? V = v f ? S3 ro 25 

i s • i , 3 ? ys-os*VM: S3 rb 25 , „„ 

Tzl; JsJ + K 1 DSsSC J)-3Cl); T F C D*DS ) 8 1 , 8 1 , 2 3 

3F=(F(.T)*F(I) ) / D S ; B = 1 ,/(OStDS + U) J U*R*DS> V=8* C 6 , *DF-V ) 

5 ’?rn=jjy fpoct) = v,* u=(2.-ij)*ds 

ys6.*3F;D3*V: TFCJ"N)2t ,31,21 

TFC< *1*2)32, 33,34 _ e 

Vs(6. *v?l-v- )/0; S3 TO 35 

V a V S J S3 TO 35 
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OSsSf J)-S(T) f IIsCpP(I)-F 0 Cn*V? FPPPCT)s(V-U)/DS 

P?P( T ) all ? FPCn = CFCJ)-FCI))/hS : DSfCV + U + U)/6. 

y=ij! J = I? I = T-K? XFC J-f041 ,5i ,41 

TsV-Kj FPPPC*i)=FPPP(n? F?PCN)=B? IN gs| 

F?C'n=DFfDi : CFPPa)d-B+B)/6.J IF(MOPE) 81 , Bi , 61 

ppopf J) =F3Mi V = F?Pf J) 

T = T ; ‘ JsJ + K i DS=SC J)-S(I) ? U*FPP( J) 

F?PPC J) =FpPP(X)+-0,5fnS*CFCn+FCJ)-DS*DS*(U+V)/12,*) 

V =U ? Iff J-M) 71,81,71 
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interoolatlon usina piecewise rAYOOR's series as generated by 
cable spline or Its integral. values f,FP,fpp,f£pp of function 
and Us derivatives are given at its nodes S fro® index JJI to 
if ihd&E > 0 a correction is added for a oiecewise constant 
fourth derivative so that integral of cubic spline is evaluated 
exactlv,i 
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5 I '3 4-1 I ' n C <4 G F*t,iSCRmf^XLSO,yLSO.XUFO,YLEL,ICM,Nt.S,MC) 
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•i L ‘4 si ST. 1M 4Llf 31 ,'VuC 31 „ _ 
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V I ^ f X ) =i\!l(.U *S0R Pf X) 1-AUC2) *XMUC 3) *X**2 

V ir> j f K l S 3 , 5 * V.J C .1 1 /SMR P ( X) *4U C 1 1 1? . * 40 C 3 ) *X 
VuTlf ) = A0C U *1DRPCXH-A[,<2>*X + 4L I ( 3m**2„ 

v t, 1 1 r ;< i = a . 5 * a u r n / s 3 r c no + a u r 2 1 + 2 , * 4 r. c 3 ) * x 
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y 4* <U = T ..,'SiJf 2 1 : Y -1 AX r . = Xl.iSuC 21 

<1 =Xf.S’U 1 1 * 
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3 2 = { fj ;I ' ! ( 7 ) - A ' H I ) * S 0 R T C X 2 1 
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x l Os* 1**7 f 
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03 10 T=1 

XCT=XCC n ; 
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FPl i Ent*0.5*A0ir/SQRTCX) + 41U 
RETURN 

♦ #*»*♦*,»!♦**♦*♦♦#♦*****•♦♦**«******♦*********#♦♦♦♦*♦♦*♦*♦***♦***♦ 
f yRCTlDtM FPLEJjl £R|OliE,XliTlil , YLTlil ,X) 


At«C3) 
, ML5TI 


C/D 

TFUCI.3T. 0.01)30 TO 102 



ill 


ADL=.-50RT 

AIL = f 
FPLELi* 
RETURN 

m 

W-tW 
*** 

fit] 


i^AOL*SQRTCXt.TLp)/XLrU 

>*AOL/SORT(X)4-Atti 





3 5 1 ,} 0 

05? 
®5? iO 

SlliS 

055 );> 

■ 5 R :)0 
06*00 
[ 0,6 ROO 
000 
’too 
! :200 
000 
'4 00 
107500 


r 

r 

c 

c 

c 

c 


10 t 


t* 4 T 3 l? 5 * 1 ?, £ hHFMF.*OP 

*1 J n R 1 ^ i f, £ .A'lFJUFCMl,!? , WJLAST»I r, A tf ) 
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SiSJS«P«3;?iSi«:gIKHiS!H!SiT«r, 

b»ck,i.j)«cikl+otkp 

TF(!C,E3,-NGF)r>V(K + l,I, JlsCIKp 

33 rn 12 

XCXS = XSIJ-XSCK-t.) 

YCYS=YGI J-YS(K-11 

RR = XFK5f*2 + YCY5tt2 , T ,, Dt . 4m ., OD 
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ADDENDUM 


The Relation Between the Number of Mean Line Elements/ Internal 
Singularities and Grid Points 


If the mean line is divided in n segments the following 
relation hold: 

Number of distributed internal sources = n 

Number of point sources (fixed) = 2 

Number of vorticity unknowns = n + 1 

Number of point vortices (fixed) = 1 

Assuming that the number of grid points chosen beyond the 

trailing edge is n , the number of grid lines in the i-directior 

w 

is 2 (n + n^) +3. if m. levels are used in the j-direction we 
have for the lifting problem/ 

Total number of field panels = Total number of grid point 

= m. x^ 2(n + r^) + 3} 

For all airfoils except the NACA 64A410 airfoil/. n w — 3. For 

NACA 64A410 airfoil/ n = 0. 

w 

For supercritical lifting flows past NACA 0012 airfoil/ 
n = 21 and for NACA 64A410 airfoil n = 22. 

For symmetric flows/ there is no vorticity distribution 

on the mean line# Hence, 

Number of distributed internal sources = n 


Number of point sources (fixed) 
Number of grid lines in i— direction 
Number of grid lines in j-direction 
Total number of field panels 


« 2 

= n + n + 2 
w 

« m 

= m x (n + n -t 
w 
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For super-critical flow past the NACA 0012 airfoil , witl 
26 x 6 grid and n = 3, we have n = 21. 

For sub-critical flows in general t the number of intern; 
singularities required are lesser since in this case close spac 
of the i = constant grid lines is not necessary. 



